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AN  ANALYSIS   OF  THREE-DIMENSIONAL   TRANSONIC   COMPRESSORS 

Antoine   Bourgeade 

February    1983 

ABSTRACT 

This  presentation  sets  forth  a  method  for  computing  the 
three-dimensional  transonic  flow  around  the  blades  of  a  compressor  or 
of  a  propeller.  The  method  is  based  on  the  use  of  the  velocity 
potential,  on  the  hypothesis  that  the  flow  is  inviscid,  irrotational 
and   isentropic . 

The  equation  of  the  potential  is  solved  in  a  transformed  space 
such  that  the  surface  of  the  blade  is  mapped  into  a  plane  where  the 
periodicity  is  implicit.  This  equation  is  in  a  nonconservative  form 
and  is  solved  with  the  help  of  a  finite  difference  method  using 
artificial  viscosity  and    artificial    time. 

A  computer  code  is  provided  and  some  sample  results  are  given  in 
order  to  demonstrate  the  influence  of  three-dimensional  effects  and  the 
blade's    rotation. 
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I.  INTRODUCTION 
Many  scientists  and  engineers  continue  to  study  the  ways  and  means 
of  making  better  use  of  the  energy  resources  at  our  disposal,  even 
though  the  energy  crisis  is  perhaps  no  longer  considered  so  severe. 
Turbines  and  compressors,  which  both  create  and  consume  energy,  have  in 
recent  years  been  the  subject  of  many  theoretical  and  experimental 
studies  aimed  at  improving  their  design  and  efficiency.  Although  some 
major  improvements  have  been  introduced  lately,  the  study  of  the 
transonic  flow  across  a  single  stage  of  a  turbine,  or  of  a  compressor, 
still  remains  highly  complex.  For  this  reason  most  of  the  theoretical 
work  done  so  far  has  focused  on  the  two-dimensional  cascade  problem 
[13-17]   or   on  the    mean    flow  problem    [22]. 

Let  us  first  discuss  the  physical  background.  A  basic  compressor 
consists  of  a  succession  of  rotors  and  stators.  These  rotors  and 
stators  are  situated  on  the  hub  and  surrounded  by  the  cowling;  they  are 
composed  of  a  certain  number  of  blades  distributed  around  the  hub,  the 
shape  of  the  blades  depending  on  their  use.  If  the  hub  is  cut  along  a 
generator  line  and  transformed  into  a  plane,  a  so-called  "cascade"  of 
these   blades    is    obtained.      A   propeller    is    a   compressor  without    cowling. 

The  present  study  deals  with  the  problem  of  transonic  flow  around 
compressor  or  propeller  blades.  From  the  mathematical  point  of  view, 
this  leads  us  to  a  system  of  partial  differential  equations  of  mixed 
type,  in  which  the  unknowns  are  the  geometric  and  physical 
characteristics  of  the  compressor.  For  a  three-dimensional  analysis 
these  equations  are  too  complex  to  be  integrated  without  some 
simplifications.  First  of  all,  we  suppose  that  the  fluid  we  are 
concerned   with   is   a   polytropic   and   nonviscous    gas   and    that      a     velocity 
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potentlal  exists.  This  necessitates  a  further  hypothesis,  namely,  that 
the  variation  of  the  entropy  is  small  so  that  the  entropy  itself 
remains  essentially  constant.  Thus  our  system  of  partial  differential 
equations  becomes  equivalent  to  a  single  equation,  the  potential 
equation,    in  both   the  steady  and    the    time-dependent    cases. 

The  potential  equation  has  already  been  solved  numerically  in  the 
three-dimensional  case  for  oblique  and  swept-wings  [2,14,15],  However, 
because  of  the  periodicity  of  the  compressor  problem,  the  square-root 
transformation  used  in  those  works  is  not  practical  here.  We  therefore 
propose  a  new  transformation,  which  maps  the  surface  of  the  blade  into 
a  plane  and    includes  periodicity   implicitly. 

The  scheme  we  use  is  similar  to  the  one  used  by  Jameson  and 
Caughey  [15]  in  the  development  of  the  swept-wing  code  known  as  FL022. 
We  solve  the  finite-difference  approximation  to  the  potential  equation 
by  row  relaxation.  A  typical  run  consists  of  50  iterations  on  a  48x6x4 
grid,  followed  by  100  iterations  on  a  96x12x8  refined  grid.  This  takes 
15  minutes  on  the  CDC  6600  and  3  minutes  on  the  STAR.  The  simplicity 
of  our  grid  generation,  together  with  the  other  hypotheses  given  above, 
limits  our  choice  of  blade  geometries.  Nevertheless,  this  method 
enables  us  to  study  how  the  speed  of  rotation  influences  the  relative 
flow  around  the  blades,  and  we  have  compared  our  three-dimensional 
results   with   data    for   two-dimensional  cascades. 

In  Section  2  we  shall  derive  the  equations  of  motion  in  physical 
space  and,  in  Section  3,  we  shall  consider  the  potential  equation, 
which  is  obtained  after  several  changes  of  variable.  The  numerical 
scheme  is  presented  in  Section  4.  The  numerical  results  obtained  are 
shown    in  Section   5   together   with  some    Calcomp   plots,      while      Section      6 
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provides      a     manual      on    the    use   of   our   computer   code.      Sections    7   and   8 
contain  a  bibliography   and   a  listing  of   the  code. 

This  study  has  been  supported  by  NASA  under  Grant  No. 
NGR-33-016-201  and  by  the  U.  S.  Department  of  Energy  under  Contract 
DE-AC-02-76ERQ3077  and  1  take  this  opportunity  of  expressing  my  sincere 
gratitude  to  Prof.  P.  Garabedlan  for  his  Invaluable  advice  and  to  Dr. 
F.  Bauer  for  her  constant  encouragement.  I  am  also  Indebted  to  my 
entire   family   for  their   moral   support. 


II.       THE    PARTIAL   DIFFERENTIAL   EQUATIONS   OF  MOTION 
This      section      sets      forth      the      various      equations      used      in      our 

computational   method    and    in   the    resulting   code. 

The      general      equations      of    fluid    dynamics    for   an    inviscid   gas   are 

well  known    [6].     These  are: 

(a)  the  equation   of    conservation   of    mass, 

(1)  pt  +  pux  +  pv+pwz+  pxii  +  PyV  +  pzw  =    0, 

where  p  is  the   density,   (u,v,w)   are   the   velocity   components,   and 
(x, y,z)  are  the  coordinates  in  physical  space; 

(b)  the  equations  of  conservation  of  momentum, 

p(ut  +  uux  +  vu   +  wuz)  +  Px  =  0 

(2)  p(vt  +  Uvx  +  vvy  +  wvz)  +  Py  =  0 

p  (wt  +  uwx  +  vw   +  wwz)  +  pz  =  ^  ' 
where  P  is  the  pressure;  and 

(c)  the  equation  of  conservation  of  energy, 


(3)  §=0  , 

at 


where   S    is   the  entropy.      It   will  be   recalled    that 


d  8  3  3.3 

=   +  u   +   V   +  w  

dt        3t  3x  3y  3z 


is     the    material    time  derivative.      We   are   neglecting  here    such  external 
actions   as    gravity. 

For  the    following   calculations    we   shall   assume   that    the   entropy   is 


-5- 

constant  throughout  the  fluid;  the  last  equation  will  therefore  not  be 
used.  But  then,  in  order  to  represent  weak  shock  waves  which  occur  in 
transonic  flows,  we  have  to  replace  the  Rankine  Hugoniot  shock 
conditions.  This  is  done  by  permitting  a  jump  in  the  horizontal 
component  of  momentum  and  by  adding  artificial  viscosity  terms  to  the 
partial  differential  equations.  The  approximation  is  adequate  for  Mach 
numbers    close   to    1. 

Let  us   suppose  now  that   the    fluid    is   a   polytropic   gas.      Therein: 

(4)  P  =    A(S)p^    , 

where  the  function  A(S),  on  the  basis  of  our  hypothesis,  becomes  a 
constant,  and  where 


is   the   adiabatic  exponent  of  the  gas,  i.e.   the  ratio  of  the  specific 
heats  at  constant  pressure  and  at  constant  volume. 

If  c  is  the  local  speed  of  sound  in  the  gas,  we  have 


(5)  c*-£. 

dp 


so  that    equations    (2)   become 

o 
p(ut  +  uux  +  vu     +  wuz)    +  c  px  =    0   , 

2 
(6)  p(vt  +  uvx  +  vv     +  wvz)    +  c  py  =    0  , 


2 
p  (wt   +  uwx  +  vw      +  wwz)    +  c  pz    =    0 
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Suppose   now   tliat    the    flow   is   steady,    i.e.       that    the    flow   is    independent 
of    t  ime ,    so  tliat 


(7)  !£.  =  o  ,    i!i=  o  ,    ll-  o  ,    £^-  o  . 

3t  3t  9t  8t 


Then  from   (6)    and    (7)   we    obtain 


(8)      p(ux  +  vy  +  wz)    -  _  (u  ux+  uvuy+  uwuz 


2  2 

+  vuv  +  v  v  +  vwv  +  wuw  +  wvw  +  w  w    )    =  0    , 


or,      if     we     assume   the   existence   of   a    velocity  potential,   <j>  ,    and    if   we 
divide   by  p/c2, 


(9)      c2(<t>xx+  «,yy+  <,zz)    -   (u2*xx+  v2*yy+  w24»zz 


+    2uv<}>xy   +    2vw<j)yz+    2wu4izx)    =   0  , 


with   u  =  <j>        v   =  <J>V,    w  =  (j>z.      This    is    the   potential      equation.         It      is 
hyperbolic   for   supersonic   flow, 

u2   +  v2   +  w2    >   c2    , 

and    elliptic    for   subsonic    flow, 

u2   +  v2   +  w2   <   c2    . 

Equations    (4)   and    (5)    show  us    that 
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(10) 


1      dcz   _    dp 


Using    the    velocity  potential,   equations    (6)   become 


1  fJ,2  Xi2,  a} 

Y 


.2 


x  +  *; +  ^)  +7rTJx  =  ° 


(ID  (y  <♦£  +  ♦£  +  ♦£)  +y^y=  °« 


so   that 


:j<*2 +  ♦£  +  ♦?)  +^rr)z  -  o  , 


(12)  I  (4>2  +  ,(,2  +  +2j    +    c      =  c0nstant 

z      x        y        z        y~i 


throughout      the      fluid.        This   is    the   Bernoulli   equation,   which  enables 
us,    knowing   the    velocity  potential,    to  compute    the    speed   of   sound. 

Suppose  next  that  the  flow  is  no  longer  steady  with  respect  to  the 
initial  frame  of  reference,  but  that,  if  we  consider  a  frame  turning 
around  the  x-axis  with  a  constant  speed  of  rotation,  co  ,  the  flow  is 
again   independent   of    time.      Then,    instead  of    (7),    we  have 


/,tN      9p  9p     ,  9p         9u  9u    ,  9u 

(15)     - —  =   -  ojz  - —  +  djy  - —  ,  - —  =   -to  z  - —  +  toy  - — 

9t  9y  9z        9t  3y  9z 


9v  9v  3v9w  9w  3w 

~ —  =   -  uz  - —  +  coy  - —   ,  - —  =   -  a)z  - —  +  wy  - — 
9t  9y  3z       9t  3y  9z 


By  introducing  the  cylindrical  coordinates  (x,9,R),  the  equations  (13) 
become 


9p  _     3p     3u  _     3u   3v  _     3v    3w  _     3w 

Ft  ~   WF3"'  Ft"  w  FT  '  Ft  ~   u  39"  '  Ft  ~  ~  w  FT 


These  new  equations  show  us  that  the  flow  dependence   on   t   and  6   is 
characterized  by  a  dependence  on  (9  -  ut)  alone. 

In  this  time-dependent  case  a  combination  of  equations  (1)  and  (6) 
gives  us 

(14)  p(ux  +  uy  +  wz)  -  2-   (u§ux+  u0vQuy+  u0w0uz  +  v0u0vx  +  v§vy 

c 

+  vQw0vz+  wQu0wx+  wQv0wy    +  w^wz)    =  0    , 

where    (uq,vq,wq)   are  the   components    of    the   relative   velocity   defined  by 

(15)  Uq  =  u    ,      vq  =    v  -  wz    ,      Wq    =  w  +  coy    . 

If    we    again  assume    the  existence   of  a    velocity  potential,    it     must      now 
satisfy   the  equation, 


(16)      c2(<frxx+t}>yy-hj,zz)-    (u^xx+  v§4>yy  +  wg<t,zz+    2u0v0<t,xy+   2v0wQ4,yz+   2w0u0<J>zx)    =   0 


This      equation      is      similar      to     equation      (9);      it      is      a   second-order 
nonlinear  partial   differential   equation  which   is   hyperbolic    for 


u2  +  v§  +  wg   >    c2    , 
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and  elliptic    for 


2  2  2  2 

15    +    VQ    +    WQ     <    c 


We  still  have  the  relation  (10)  for  the  speed  of  sound,   so   that 
using   the  velocity   potential   and   with  the  help  of  (13),  we  obtain, 
instead  of  (6) , 

1  fj.  2^  2^  2 


(17)        (_  ((fr^-h^)  -a)Z4»y+o)y*z  +^IT)X  =  0  , 


and  this    gives  us    the   new  Bernoulli   equation, 


(1«)       y  ^x^y4*^    "  wz<f,y  +  ^z  +  4rf  =  Constant 


along   each    line    parallel    to   the    x-axis. 

In    the   next   section  we      shall     transform     these      equations,      by      a 
change   of  coordinates,    into    the    system  which   is   solved  by  our   code. 
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III.      GRID  GENERATION 

We  have  obtained  the  equations  to  be  solved,  but  we  have  still  to 
impose  the  geometrical  constraints  due  to  the  slip  condition  on  the 
surface  of  the  compressor  blade  and  to  the  periodicity  condition.  In 
this  section,  we  shall  describe  the  geometrical  space  in  which  the 
equations  will  be  solved.  We  shall  therefore  list  all  the  mappings 
which  are  performed  to  transform  the  physical  space  onto  the 
computational    space. 

For  the  purpose  of  the  periodicity  condition,  we  begin  by 
introducing  the  angle  9  of  the  cylindrical  coordinates  (x,G ,R)  with 
respect    to    the   x-axis.      To    accomplish    this    we  use   a    conformal  mapping 

(19)  (x,y,z)  -    (xo,0 ,ZQ) 
defined   by 

(z+iy)    =  exp(Z0+i8),      x  =  xn    . 

The  Jacobian  matrix  associated  with  this  transformation  is 
1   0   0 

(20)  Jl   =  0   P   Q 

0     -Q        P 
with   elements    defined  by 


_        cos  6  sin  9  ,„    v 

P  =    — —  ,      Q  =  — - —  ,      R  =  exp(Zn)    . 


We    also   require    the    inverse   matrix,   namely, 
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1   0   0 
JJ1  =       0   z  -y 

0    y    z 

->■->■  + 
If   (i,j,k)  is  the  orthonormal  basis  for  physical  space,  then  the  basis 

->■   +■      ■* 
vectors  connected  with  the  new  coordinates  are  i,  u,  and  v,  where 


(21)  u  =  zj  -  yk  ,     v  =  yj  +  zk  . 

For  the  surface  of  the  blade,  the  representation  of  the  finite 
difference  slip  condition  becomes  greatly  simplified  and  more  accurate 
if  the  boundary  surface  lies  on  a  coordinate  plane.  The  idea,  cf. 
[14],  is  to  transform  the  surface  of  the  blade  onto  a  plane  which  will 
constitute   the    lower  boundary   of  a   half-space. 

Unfortunately,  if  we  apply  the  square-root  transformation  of 
reference  [14] ,  the  periodic  strip  is  transformed  in  such  a  way  that 
the  periodicity  condition  is  difficult  to  satisfy.  We  therefore  prefer 
a  transformation  which  would  map  a  periodic  strip  conformally  onto  a 
half-space,    so  that    the   periodicity    condition   becomes   implicit. 

The  required  transformation  can  be  decomposed  into  two  successive 
mappings.  For  the  purpose  of  simplification  we  consider  only  plane 
sections  orthogonal  to  the  ZQ-axis.  The  first  mapping  transforms  a 
periodic  strip  (cf.  Figure  la)  onto  a  slit  plane  (cf.  Figure  lb). 
The  image  of  the  two  lines  delimiting  the  strip  is  the  negative  real 
axis.  We  then  apply  a  square-root  transformation  which  maps  this  plane 
onto   a  half-plane   (cf.      Figure   lc). 
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Figure  la.   The  periodic  strip. 
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*      C 


Figure  lb.   The  slit  plane, 


<£ 
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Figure  lc .   The  half  plane  with  a  bump. 
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To    Implement    these    transformations   we  draw  in  each  plane    section  a 
singular      line      (the     hranch   cut)   from   inside   the   blade,    near    the   nose, 
out    to   downstream  infinity.      This   is    actually  a  half-line  whose      origin 
has    the   coordinates    (xg,6s)    .      We    next    perform   the  mapping 

(22)  (x0,9,Z0)  +    (X.Y.Z)    , 
defined  by 

2 

(x0-xs)  +  i(e-es)  -n"1  logfi  +  (X+*Y)  )   ,    z0  -z  , 

where  N  and  T  are  two  constants.  N  is  more  precisely  the  number  of 
blades  on  the  compressor.  The  coordinates  xg  and  6g  usually  depend  on 
Z.  Their  derivatives  with  respect  to  Z  will  be  denoted  xz  and  9Z  , 
respectively. 

The  Jacohian   matrix  determined  by  this    transformation   is 
a     -b        0 

(23)  J2   =     b        a        0 

e        f        1 
whose  elements    are    defined  by 


2  2 

a    =  Hxx    ,    b    =  H6X   ,    e  =    -  axz  -  b6z    ,      f  =    -  a6z  +  bxz    ,      H  =  xx  +  9X 


The  inverse    matrix   is 
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XX        8X  ° 

J  2      =  ~°X        XX  n  ' 

xz     e  z       i 

and    the    basis    vectors    related    to    this    transformation   are      now      (4,B,C), 
with 

(2  4)      A   =   Xyi   +  9vii       ,       B  =    -  O^i   +  x^u    ,       C  =   x?  i   +  6  7u   +   v 

Let         us         return      to      the      potential      equation      (9),      it      can      be 
reformulated   as 

(2  5)  c2V24,    -   (V4>    •    V)2<f>    =   0   , 


where 


V    =  —  i  +  —  j   +  —  k 
9x  d y  3z 


is   a    notation    for   the    gradient.      If    we   denote    the    final   Jacobian   matrix 


A  =  J  i    •    J2 


hv 


and,    if   we  use 


A     =      (a±J) 


9  A  9  A        -      d 
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for  the    derivatives   with    respect      to      the      new     coordinates,      then      the 
derivatives   with    respect    to   the   initial  coordinates   are    represented  hy 


Jt=\   aUdJ    '      37  =  £    a2jd.i    •  J^  =  l    a3jdj 
J  J  j 

Thus  the  Laplacian,   which  is    defined  by 


v2  =li  +  li+  li 

9x2        3y2        3z2 


can  be   rewritten  as 


(26)      V2    -       I        aijaikd.jdk  +       £        aijaj!kdk.    » 


where      the      coefficients      a£  .  stand    for   the    derivatives    of    the  Jacobian 
matrix  elements; 


lik   =  dj  aik 


Let  us   set 


(Fi.)    =  F   =   AC«A    ,      DL    =  V2X    ,       D2    =  V2Y    ,      D3    =  V  2Z 


These   notations   allow   us   to   simplify    the   equation   (26)    to  arrive   at 


(27)  V2    =     I      Fj    k   djdk  +  I    Dkdk    . 

j,k  k 


For  the  second  term  of  equation  (25)  let  us  denote  the   components 
of  the  velocity  in  the  physical  space  by 
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9<}>  di>  3$ 


Using   the    same   computation  as    for   the   Laplacian,   we  obtain 

(28)  (V^-V)2    =     I       \B£dkdJl   +  I    C^d^    , 

k,£  I 

->•-»•->■ 
where   (B^Bp.B^),  the  components  of  the  velocity  in  the  basis  (A,B,C) 

are  given  by 


B.  =  I    a1J+1 


and  where  the  coefficients  C.    are  defined  by 


C1  =  (V-j).V)2X  ,   C2  =  (V^j.V)2  Y,   C3  =  (V^V)2  Z 


After  these  necessary  but   somewhat   tedious   transformations   the 
potential  equation  acquires  the  useful  form 

(29)     I       (c2Fti  -  BiB^d^jO.  +  I    (c\    -  C±)    dtf   -  0  . 
i.j  1 

This  has  the  advantage  of  being  relatively  tractable,  for  an  equation 
which  is,  really,  quite  complicated.  In  the  time-dependent  case,  we 
need  only  to  replace  the  cf>.'s  by  the  components  of  the  relative 
velocity  in  the  calculation  of  the  coefficients  of  equation  (29). 
Moreover,  if  a  reduced  potential  is  defined  as  the  difference  between 
the  true  potential  and  the  uniform  flow  potential  corresponding  to   the 
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inlet      speed,    we  can   use    it    to  compute   the    derivatives    d . d  «4>    and  d.<|>    in 
equation  (29)  without    changing   that   equation. 

We  next  describe  the  grid  used  in  our  computational  method.  It  is 
obtained  with  the  help  of  a  system  of  sheared  coordinates.  These  are 
defined  by  considering  coordinates  parallel  to  the  transformed  surface 
of    the   blade.      If 

Y  =    S(X,Z) 

is      the      equation      of      this      surface,    the    transformation   in   question   is 
defined   by  setting 

T   =   X    ,      J  =    Y   -    S(X,Z)    ,         K  =    Z    . 
The  Jacobian   matrix   is 

1      -sx        0 

(30)  J3  =  0  1  0 

0      -sz        1 

and  the  related  basis  vectors  are 


(31)    I  =  A  +  SYB  ,    J  =  B  ,   K  =  C  +  S7B 


In  order  to  carry  out  the  computations  on  a  finite  domain,  these 
last  coordinates  are  stretched  so  that  the  final  domain  of  computation 
becomes  a  cube  the  edges  of  which  have  length  1.  After  these 
transformations  equation  (29)  remains  of  the  same  general  type,  but  its 
coefficients  have  to  be  changed  slightly.  Complete  details  are 
specified  in  the  listing  of  the  code  in  Section  VIII. 
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It  remains  to  define  the  boundary  conditions  at  the  hub,  at  the 
cowling,  and  on  the  blade  surface.  We  impose  a  slip  condition  for  the 
flow  on  these  boundaries;  this  leads  to  a  Neumann  problem  for  the 
velocity  potential,  since  the  normal  derivatives  at  the  boundaries  must 
vanish.      We   shall   express    this    condition    in    our    system   of    coordinates. 

The  hub  and    the    cowling   are    defined  by  equations    of    the    form 

Z    =  Constant . 

■*■-*■  > 

At    any    boundary   points    the    two    tangential    vectors    are    i    and   u,    and    v 

>        -* 

i    x    u    is    a    normal    vector.       The    boundary    condition   at    the    hub  and    at    the 

cowling   may   thus    be  expressed  by  the   equation 


(32) 


3<J> 


y<J>  9    +  Z(| 


On    the    other    hand ,    the    blade   surface   is    defined  by   the   equation 


J   =  0 


The   two    tangent    vectors    I   and  K    lead    to    the    normal    vector   I    x       K     whoso 


coordinates    in   the   basis   (i,u,v)   are 


a   =  R"(9K  +  sxxx)    '       8  =   SX8X"XX    '   y    = 


Sve  -  f  +  S7 


Hence  the    boundary   condition  on    the    blade   surface    is    given  by 
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(33)      ex  P-  +  &  |1  +Y   #--«♦!   +  ©*  +  YYH2   +   (Tz-6y)«l»3   =   0 

0  Xa  0  0  O  £rv 


In  equations  (32)  and  (3  3),  (<j>  ^  ,<j>  ^ » 4»  3)  must  be  the  derivatives  of 
the  true  potential,  i.e.  the  components  of  the  velocity.  For  the 
time-dependent  case,  they  are  replaced,  in  equation  (33),  by  the 
components  of  the  relative  velocity. 

In  the  code  the  values  of  (ij>  ■  ,<j>  jA  3)  are  computed  only  when 
necessary,  and  the  only  derivatives  available  are 

u  =  i*      v=il,   w  -  li  . 

ax  '  ay  '  az 

These   derivatives  are   related    to   the   components    of    the    velocity   through 
the   relation 

(34)  4>2        =   A        V 

*3  W 

The  boundary   conditions    (3  2)    and    (33)    become 

(3  5)  eU  +   fV   +  W  =   0 

and 

(36)      (oa  +  6b  +  ye)U  +   (-ob  +  Ba  +  yf)V   +  yW  =    0  . 

In  three  dimensional  space  an  appropriate  vortex  sheet  behind  a 
blade  must  be  considered.  The  shape  of  the  vortex  sheet  is  modeled  by 
our  code  in  the  following  way.  In  each  plane  cross-section,  the  upper 
and   lower  surfaces  of  the  blade  are  extended  behind  the  trailing  edge 
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by  two  lines  parallel  to  a  branch  cut  (Cf.  Figures  1,4,7).  These  two 
lines  represent  the  vortex  sheet  in  that  section,  and  if  the  blade  is 
closed  they  are  identical.  For  two  points  situated  on  opposite  sides 
of  the  vortex  sheet  the  pressures  are  the  same,  and  the  normal 
velocities  are  zero;  only  the  tangential  components  of  the  velocity  may 
be  different.  Since  a  shape  is  assumed  for  the  vortex  sheet,  we 
require  only  continuity  of  the  normal  component  of  the  velocity  across 
the   sheet.      Computationally  this    reduces   to  a    condition    like 

(37)  (frYY  =  0 

after  a  jump  in  <j>  is  removed.  Moreover,  the  jump  of  the  potential  $ 
across  the  vortex  sheet  is  supposed  to  be  constant  in  each  plane 
section.  The  Kutta-Joukowsky  condition  at  the  trailing  edge  determines 
this    jump.      This  amounts    to  a    linearized    treatment    of   the    vortex   sheet. 
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IV.       FINITE    DIFFERENCE    SCHEME 

The  success  of  codes  for  the  design  and  analysis  of  supercritical 
wings  [1,2,3,15]  shows  how  effective  the  computational  fluid  dynamics 
has  become  for  transonic  flow.  The  first  step  in  this  development  was 
the  design  of  shockless  airfoils  by  the  hodograph  method  [1].  Then  the 
introduction  of  a  retarded  difference  scheme  [20]  allowed  the  analysis 
of  flow  at  off-design  conditions.  This  scheme  incorporated  artificial 
viscosity  in  order  to  capture  shocks  in  the  supersonic  zone.  It  was 
then  improved  to  permit  the  analysis  of  three-dimensional  wings  [14] 
and  their  design  [11].  In  our  computational  method  we  use  the 
last-mentioned  scheme,  and  this  section  explains  how  it  is  incorporated 
into   our   computer   code. 

All  the  equations  we  have  described  in  the  previous  section  are 
represented  in  the  computer  code  by  finite  difference  approximations. 
The  flow  conditions  at  each  grid  point  are  determined  with  the  help  of 
the  reduced  potential  C.  The  first  derivatives  of  this  potential  are 
calculated  by  using  central  differences.  With  these  values  we  compute 
the  approximate  velocity.  This  allows  us  to  determine,  with  the  help 
of  Bernoulli's  equation  (12),  whether,  at  the  point  considered,  the 
flow    is   subsonic  or    supersonic. 

If  the  flow  is  subsonic  the  second  derivatives  used  for  the 
computation  of  the  potential  equation  residual  are  approximated  by 
central   differences    of   the    form 


(38) 


Gi+l,j,k        2Gi,j,k  +  Gi-l)j,k 
(AX)2 
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and 


(3  9)         0 


'XY 


Gi+1,  j+l,k  "  Gi-H,j-l,k  "  Gi-1,  j+l,k  +  Gt-1,  j-l,k 
_______ 


For  supersonic  points,  on  the  other  hand,  we  have  to  introduce 
artificial  viscosity  in  order  to  capture  weak  shocks.  This  is 
accomplished  by  using  a  retarded  difference  scheme  [15].  Thus  we 
separate    the   equation    (29)    into    two  groups    of    terms: 


(40) 


U2-q2)  Gc.  +  c2(V2G     -  G„) 


where  q  is  the  speed  and  s  a  coordinate  in  the  flow  direction.  The 
first  term  represents  the  second  derivatives  in  the  flow  direction;  the 
derivatives  in  the  other  directions  form  the  second  term.  The  second 
derivatives  used  in  this  second  expression  are  computed  by  using 
equations  (38)  and  (39).  But  for  the  second  derivatives  of  the  first 
expression  we  use  retarded  difference  approximations  of  the  following 
types: 


(AX)' 


a  nd 


(42)  G 


rovided    that    the    velocity   has    positive   components    in    the      X      direction 
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and  in   the   Y    direction.      Equations    (41)    and    (42)    are      only     first-order 
accurate   and    introduce  the    truncation  errors 


-  AVG 


AXGyvy   "*"  AY    GyvY 


XUXXX   ' 


For  the  potential  equation  (40),  at  supersonic  points,  these  terras 
represent  a  positive  artificial  viscosity  which,  when  the  flow  is 
aligned  with   the   X  direction,    reduces    to 


(q2    -   c2)   AX    G^x 


as    in    the    scheme   of   Murman  and   Cole    [20], 

Equations  (38)  to  (42)  are  used  at  each  iteration  of  a  run.  To 
describe  the  iteration  process,  which  is  done  by  row  relaxation,  it  is 
helpful  to  introduce  an  artificial  time  t,  which  increases  by  the 
quantity  At  at  each  iteration.  The  right-hand  sides  of  equations  (38), 
(39)   and    (41)    become 


G?+l,j,k-    ^-^   Gij,k  -    (2/">   Gi,j,k+  Gi-1,  j,k 
(AX)2 


Gi+l,j+l,k"  Gi-l,j+l,k  "Gi+ltj-l)k+  Gi-1, j-l.k 

4AX  AY 


2Gi,j,k-  Gi,j,k-    2Gi-l,j,k+  Gi-2,j,k 
(AX)2 


where      the      superscripts      N      and     0      denote      new     or     old   values    of   the 
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potenttal      and    where  w     is    the   overrelaxat ion    factor.      These   expressions 
represent  approximations    to 


GYY  -  —  (GYt.    +  —   (-  -    DG.),    GYY   -  I—  GYt    ,      GYY  +   2  ^  G, 
xx       AX        Xt        AX     (o  fc    •       XY        2  AX     Yt  xx  AX      > 


Hence   the   equation    solved    is 

(43)        I      (cVr  B^GX^   +  ?    1<CV  C0%  +  aiGXit]    +6Gt   -    °   • 
i,j  l 

where      the      coefficients     a,      and      5      are        determined        by        the        new 
approximations    (38),    (39)    and    (41). 

By    considering   an  orthonormal    system  of    coordinates    (p,r,s),    where 
the    flow  direction   is   still   the    s    direction,   we   arrive   at    the   equation 

(4)       (c2-q2)Gsg   +  q2Grr+q2Gpp    +6^+  B2Grt+  ^pt  +  5  Gt    =    °   ' 
If 


1  r     Pi  $2  ^3     , 

t    =   t    -  _     s   +  —  r  +  —  p 

-      c    -q-  q"  q" 


is  a  new  time  coordinate,  then  equation  (44)  is  transformed  into 


2     n  ?    o  2 
(45)   (c2-q2)Gss  +  q2Grr  +  q2G   -  (_— -  +  —  +  —  )GTT  +  5GT  =  0. 

c  -q    q    q" 


In    order   to   ensure      the      convergence      of      the      iteration      scheme      to      a 
solution     of      the      steady-state  equation,   we  want    equation   (45)    to  be   a 
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daraped    three-dimensional    wave    equation.      Given      Che      form      of      equation 
(45)    this    is    equivalent    to   the   condition 


2  2  2 

(46)  (-T2  2       -T)(c      q   >    >  °  ' 

c    -q-         q  q 


At  subsonic  points  the  damping  condition  is  always  satisfied.  At 
supersonic  points  it  may  no  longer  be  true.  The  choice  of 
approximations  (41)  and  (42)  help  to  ensure  a  large  coefficient  g,  , 
but  near  the  sonic  line  this  may  not  be  enough.  One  way  to  ensure  that 
condition  (46)  is  satisfied,  is  to  increase  g,  by  adding  a  term  of  the 
form 


6  U  °st  "  e  H  ((flGxt+  *2Gyt  +  *3Gzt>  =  ^  ||  (*iGXt+  B2GVt+  B3GZt) 


to  equation  (43),  where  (3  is   a   positive   parameter   large   enough   to 
increase  B  j  by  the  right  quantity. 

The  mixed  derivatives  used  above,  G„   ,  Gy,-  ,  G7r  are  obtained  by 
using  approximations  of  the  type 


rN       P0       PN       ,  p0 
,   =  Gi,j,k  -  Gi,j,k  ~  Gi-l,j,k  +  Gi-l,j,k 

'Xt  At  AX 


We  have  described  the  difference  scheme  for  interior  points.  For 
the  boundary  conditions  at  the  hub  and  cowling,  we  transform  equation 
(35)  into  a  difference  equation  in  order  to  compute  the  value  of  the 
reduced     potential      at  these   boundaries.      For  the   blade   surface  and    the 
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vortex  sheet  we  introduce  ghost  points  behind  the  boundary.  These 
points  are  used  to  compute  the  value  of  the  potential  at  points  on  the 
blade  as  if  they  were  interior  points.  The  values  of  the  potential  at 
the  ghost  points  behind  tlie  blade  are  determined  by  using  equation 
(36),  i.e.  the  slip  condition.  For  the  points  behind  the  vortex  sheet 
the  potential  is  determined  by  the  Kut ta-Joukowski  condition.  To 
compute  the  value  of  G  at  the  remaining  boundary  points  we  use  an 
approximation    to    the    outlet    velocity. 

The  program  permits  not  only  compressor  blades  but  also  propellers 
to  be  modeled.  Since  propellers  have  no  cowling,  the  inlet  speed  is 
equal  to  the  outlet  speed.  Thus,  for  a  propeller  run,  all  the  boundary 
conditions  except  those  at  the  hub,  on  the  blade  and  on  the  vortex 
sheet   are    replaced   by 

G   =  0    . 

With  program  FL022  as  the  starting-point,  we  have  used  the  result 
stated  in  this  section  to  write  a  new  program  in  FORTRAN  named  CSCDF22. 
This    program  is    listed    in   Section   Vltl. 


-28- 

V.      RESULTS 

In  the  preceding  sections  we  described  our  computational  method. 
In  order  to  show  how  tliis  method  works,  we  present  in  this  section 
seven  examples    of    runs   made   with  the   program  CSCDF22. 

The  first  case  was  run  to  test  the  validity  of  our  code.  A 
compressor  blade  had  been  designed  with  code  K  [3].  Figure  2  shows  the 
result  of  this  run.  We  used  the  coordinates  of  this  two-dimensional 
blade  profile  to  create  a  three-dimensional  blade.  The  profiles  were 
identical  at  each  span  station.  Since  we  chose  to  run  the  program  for 
a  cascade  configuration  and  in  the  compressor  case,  we  were  solving  a 
two-dimensional  cascade  problem.  The  inlet  speed  and  the  inlet  angle 
were  the  same  as  in  the  design  run,  but  the  blade  was  markedly  cambered 
and  we  chose  a  gap-to-chord  ratio  of  1.5.  Figures  3  and  4  show  a 
representation  of  the  cascade  and  the  grid  at  the  hub  which  is  the  same 
at   each   span   station.      The  pressure   distribution   is    given   in  Figure   5. 

The  other  six  examples  are  all  run  with  the  same  blade.  The 
profile  is  defined  at  three  different  span  stations.  The  basic  profile 
is  the  NACA-0012  profile.  At  each  span  station  the  chord  and  thickness 
are  the  same  but,  for  the  realistic  axial  flow  configuration,  the 
coordinates  of  the  singular  points  are  different.  The  blade  has  a 
sweep  angle  of  14  degrees  and  the  dihedral  angle  decreases  from  0. 
degrees  at  the  hub  to  -10.  degrees  at  the  tip.  The  angle  of  twist  is 
equal  to  0.  degrees  at  the  hub,  2.  degrees  at  the  cowling,  and  1. 
degree  in  between.  For  the  axial  flow  configuration  the  distance 
between  the  hub  and  the  axis  is  2.  This  is  also  the  distance  between 
the  hub   and    the    cowling    (or  the    tip). 

The  first   example  with   this      blade      is      a      run      for     a      compressor 
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cascade  configuration  which  appears  in  Figure  6.  Figures  7  and  8  show  a 
representation  of  the  grid  and  the  pressure  distribution.  The 
different  flow  parameters  are  shown  in  Figure  8.  Ml  is  the  inlet  Mach 
number.  M2  is  the  outline  Mach  number.  DEV  is  the  difference  between 
the   outlet  angle   and    the    inlet   angle    given  by  ALP. 

The  four  ne Kt  runs  are  also  analyses  of  compressor  flow,  but  now 
in  an  axial  flow  configuration  and  with  eight  blades  around  the  hub. 
For  each  run  the  parameter  OM  determining  the  speed  of  rotation  has 
different  values.  These  values  are  0.0,  0.5,  1.0  and  2.0.  The  other 
flow  parameters  are  identical.  Figures  9  and  10  represent  the  grid  at 
the  hub  and  cowling.  It  will  be  noticed  how  different  they  are, 
although  the  profile  is  almost  the  same  at  the  hub  and  cowling;  but  the 
gap-to-chord  ratio  is  twice  as  large  at  the  cowling.  Figure  11  shows 
how  the  blade  looks  in  a  plane  orthogonal  to  the  axis.  The  different 
pressure  distributions  obtained  for  each  run  are  given  in  Figure  12, 
13,  14  and  15.  These  figures  clearly  illustrate  how  the  flow  evolves 
as  the  speed  of  rotation  increases:  the  shock  on  the  upper  surface 
weakens  and  then  disappears  while,  on  the  other  hand,  a  shock  appears 
on    the    lower   surface  and    intensifies. 

The  last  example  is  similar  to  the  fifth,  except  that  it  is  a 
propeller  analysis.  Figure  16,  17,  and  L  8  show  the  geometry  of  the 
case,  and  the  pressure  distribution  is  given  in  Figure  19.  by 
comparing  Figures  14  and  19,  we  observe  that,  for  example  5,  the  shock 
appearing  at  the  cowling  which  is  caused  by  the  speed  of  rotation,  is 
amplified  significantly   by   the    cowling. 

These  numerical  experiments  show  that  our  computat ional  method,  in 
spite      of      its      restrictive      hypotheses,      enables      us      to      analyze      the 
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th  ree-dimens  tonal  transonic  flow  around  compressor  blades.  This  Is  a 
first  step  towards  the  study  of  compressors  or  propellers  in  three 
dimensions.  Considering  what  happened  in  the  study  of  transonic  flow 
past  swept  wings,  we  venture  to  suggest  that  the  next  steps  could  be  to 
improve  this  method  in  order  to  include  the  computation  of  the  wave 
drag  and  to  design  blades  with  a  prescribed  pressure  distribution 
[4,11]. 
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1.2  1  CP 
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MU.707         M2=.534       DEL  TH=  35-00       G/Cr    .99 


Figure   2 


-32- 


CflSCRDE    REPRESENTATION 


G/C    =    1.  50 


Figure  3 
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G-RID  ON  THE  SURFflCE   Z  =  0.  00 


Figure  4 
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PRESSURE    DISTRIBUTION  BLADE    PROFILE 

Ml    =      .71,        M2    =         .57,        DEV    =    -32.0.       ALP    =    35.0 


Figure  5 
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CASCADE   REPRESENTATION 


G/C   =    1.50 


Figure   6 
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GRID   ON  THE   SURFACE      Z   =   0.00 


Figure  7 
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PRESSURE   DISTRIBUTION  BLADE   PROFILE 

M1    =      .74,      M2   =        .70,      DEV  =  .2,      ALP   =      0.0 


Figure   8 
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GRID    ON    THE    SURFRCE      Z    =    2. 00 


Figure  9 
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GRID  ON  THE  SURFACE   Z  =  4. 00 


Figure  10 
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CUT    IN    THE    PLRNE    X    =    . 500 


Figure  11 
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PRESSURE    DISTRIBUTION  BLADE    PROFILE 

Ml  =   .75.   M2  =    .  71.   DEV  -     .2,   ALP  =   0.  0 


Figure  12 
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PRESSURE  DISTRIBUTION      BLADE  PROFILE 
Ml  =   ■  75,   M2  =    .  70,   DEV  =    1.  9,   RLP  =   0.  0 


Figure  13 
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PRESSURE    DISTRIBUTION  BLADE    PROFILE 

Ml    =      .75.       M2    =         .70.       DEV    =         3.5.       ALP    =      0.0 


Figure  14 
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PRESSURE    DISTRIBUTION  BLADE    PROFILE 

Ml    =      .75,       M2    =         .69.       DEV    =         7.0.       RLP    =      0.0 


Figure  15 
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GRID  ON  THE  SURFRCE   Z  =  2. 00 


Figure  16 
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&RID    ON    THE    SURFnCE      Z    =    4. 00 


Figure  17 


-47- 


CUT    IN    THE    PLANE    X    =    . 500 


Figure  18 
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PRESSURE    DISTRIBUTION  BLADE    PROFILE 

Ml    =      .75,       M2    =         .75.       DEV    =         O.O,       ALP    =      0.0 


Figure  19 
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VI.       HOW   TO   USE   THE   CODE 

This  section  is  intended  to  serve  as  a  guide  for  users  of  our 
computer  code.  We  explain  the  different  code  options  and  give  a 
listing   of    the    input    parameters. 

Our  computer  code  includes  two  options  which  permit  the  solution 
of  either  the  compressor  problem  or  the  propeller  problem,  and  it 
treats  them  either  in  a  realistic  axial  flow  configuration  or  in  a 
three-dimensional  cascade  representation.  In  the  first  option  there  is 
a  choice  of  two  boundary  conditions.  In  the  propeller  case  the  reduced 
velocity  potential  is  set  to  zero  at  infinity;  in  the  compressor  case 
the  boundary  condition  is  given  by  a  Neumann  condition  at  the  cowling. 
The  first  option  determines  the  first  mapping  to  be  performed.  For  the 
axial  flow  configuration  the  whole  space  is  mapped  into  a  part  of  the 
space  consisting  of  as  many  periodic  strips  as  there  are  blades  around 
the  hub.  For  the  cascade  configuration  this  first  mapping  is  not 
required,  since  the  computation  starts  directly  from  a  periodic  strip. 
In  both  cases  the  periodic  strip  is  mapped  conformally  into  a 
quadrilateral.  For  the  compressor  case  the  blade  shape  and  the  vortex 
sheet  define  the  bottom  of  the  quadrilateral;  the  hub  and  the  cowling 
determine  the  sides;  and  the  upper  edge  is  defined  by  the  flow  at 
infinity.  For  the  propeller  case  the  side  previously  determined  by  the 
cowling   is  defined  by  the    flow  at  infinity. 

The  shape  of  the  blade  is  defined  by  giving  the  coordinates  at 
different  span  sections  from  the  hub  to  the  cowling.  As  many  as  six 
sections  can  be  defined  by  giving  the  Cartesian  coordinates  at  each 
section  with  z  constant.  If  two  sections  are  similar,  only  the 
coordinates    of   the    first   section  are   read   in. 
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We    shall   now     describe      the      input      data      cards      and      explain      the 

parameters      which      occur      in   the   input.      All  the    data    cards  are    read   in 

format   (8F10. 6).      A  title  card    is   inserted  before      each      data      card      or 

each   series    of    data    cards. 
Title    Card   1: 
NX,   NY,    NZ,    FPLOT,    XSCAL,    PSCAL,    GRIDX,    GRIDY. 

NX,    NY,    NZ  are      the      number      of   mesh   cells   in  each   direction  of   the 

transformed  space.  Thus  NZ  is  the  number  of  mesh  cells 
in  the  radial  direction,  NX  the  number  of  mesh  cells 
along  the  blade  and  the  vortex  sheet ,  and  NY  the  number 
of  mesh  cells  in  the  third  direction  of  the 
computational  grid.  If  NX  =  0 ,  the  program  stops.  The 
dimensions    in   the   code  are  NX+1  ,   NY+2,   NZ+3. 

FPLOT  is     the      parameter      controlling    the    generation  of   plots. 

IF  FPLOT  =  0,  a  printer  plot  but  no  Calcomp  plot  is 
obtained  at  each  span  station.  If  FPLOT  =  1,  both  a 
printer  plot  and  a  Calcomp  plot  are  generated;  and  if 
FPLOT  =   2,   only  a    Calcomp   plot   is    generated. 

XSCAL,  PSCAL  determine  the  scales  of  the  Calcomp  plot.  The  pressure 
scale  is  set  to  PSCAL  per  inch  in  each  section  plot. 
PSCAL  =  0  is  equivalent  to  PSCAL  =  0.5.  If  PSCAL  is 
positive,  each  section  is  scaled  so  that  the  span  is  5. 
If  PSCAL  is  negative,  each  section  is  scaled 
proportionately  to  the  local  chord  and  the  maximum  chord 
is  XSCAL/2. 

GRIDX,  GRIDY  control  the  grid  generation.  The  grid  depends  largely 
on     the      geometry     of      the      compressor,      so      these        two 
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parameters  must  be  chosen  carefully  so  as  to  obtain  an 
acceptable  grid.  GRIDX  has  to  be  greater  than  -1  and 
GRIDY  is  positive.  In  most  of  the  cases  we  can  use 
GRIDX  =  0  and  GRIDY  =  1.  Decreasing  GRIDX  gives  more 
points  near  the  nose  and  less  at  the  trailing  edge. 
Increasing  GRIDX  gives  fewer  points  near  the  nose  and 
more  at  the  trailing  edge.  Increasing  GRIDY  increases 
the  number  of  points  near  the  blade  and  the  vortex 
sheet . 

Title   Card   2: 

MIT,    COV,    PI,    P2,    P3,    BETA,    FHALF ,    FCONT. 

MIT  is  the   maximum  number  of   iterations    computed. 

COV  is    the    desired   accuracy.      If    the   maximum      correction      is 

less      than     COV,      the      program     goes    to   the  next   grid  or 

terminates.        The      same      will     happen      if      the        maximum 

correction  is   greater   than   10. 
P1,P2  are      the      relaxation   parameters   for   the    subsonic   and   the 

supersonic   points    respectively   PI   is  between    1.    and      2., 

P2  is  less    than  or   equal    to   1. 
P3  is     used      for   the   boundary  values.      P3  is   usually  set   to 

1,    but  a    smaller   value   improves    the   convergence    for  very 

cambered  blades. 

BETA  is  the   damping   parameter   controlling    the  amount    of   added 

d>    ...      It    is  normally  set    between  0.    and  0.5. 
T  s  t  J 

FHALF  indicates   whether    the   mesh  is    to  be   refined.      For     FHALF 

=     0      the     program     stops   after   convergence   or   after   the 
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FCONT 


prescribed  number  of  iterations  have  been  completed. 
For  FHALF  =  1  the  mesh  will  be  refined  and  another  input 
card  must  be  read  for  the  parameters  of  the  refined 
mesh.  The  mesh  can  be  refined  only  twice  and  the 
maximum  grid   is    128x  24x14. 

indicates  how  the  computation  starts.  After  each  run 
the  velocity  potential  is  stored  on  tape  8.  If  FCONT  =  0 
the  program  begins  by  initializing  the  velocity 
potential.  If  FOONT  =  1  ,  we  have  a  continuation  run  and 
the  data  stored  on  tape  8  from  a  previous  run  is  used  as 
input  and   is  read  in  on   tape   7. 


Title   Card  3: 

FMACH,    OM,    AL,    PA,    ZONE,    DC. 


FMACH 
CM 


AL 
PA 


ZONE 


is  the   inlet   Mach  number. 

is        the        speed      of      rotation.        This     means      that      the 

compressor   is   turning  at  the   speed  of   OM  revolutions  per 

60     units    of   time.      The  unit   of   time  depends  on  the   flow 

conditions.      The  unit   of    length  is      the      length     of      the 

blade      divided      by     CHORD     and      the   unit    of   speed    is   the 

inlet   speed. 

is  the  inlet  angle. 

determines  whether  an   axial   flow  configuration  or   a 

cascade  is  to  be  studied.   PA  =  0  for  a  cascade.   PA  =  1 

for  an  axial  flow  configuration. 

is  the  distance  between  the  axis  and  the   hub.    In  the 

cascade  case  it  is  set  equal  to  1. 
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DC  is      the    gap-to-chord  ratio    for   the   cascade   case.      For  an 

axial  flow  configuration  it  is  the  number  of  blades 
around    the   hub. 

Title   Card   4: 

ZSYM,    NC,    SWEEP1,    SWEEP2,    SWEEP,    DIHED1,    DIHED2,   DIHED 

ZSYM  indicates      whether     we     have     a      compressor   problem  or  a 

propeller  problem.  ZSYM  =  0  for  the  compressor  case  and 
ZSYM  =  1  for  the  propeller  case.  Thus,  if  ZSYM  =  1, 
there  is  no  cowling,  the  outlet  speed  is  the  same  as  the 
inlet  speed  and  the  radial  direction  is  stretched  to 
infinity. 

NC  is  the  number      of      span      stations      at     which      the      blade 

section  is  defined.  The  blade  section  is  interpolated 
linearly  between  two  span  stations.  If  NC  <  3  the 
geometry  is  given  by  the  last  case  studied  and  the 
program  will  run  for  this  ZSYM  and  for  the  parameters 
given  by  data  cards  1  through  3  without  the  need  to  read 
any  more   cards. 

SWEEP1 .SWEEP 2, SWEEP  are  the  angles  of  sweep  of  the  singular  line  at  the 
hub,  at  the  cowling  (or  the  tip)  and  at  the  far  field 
respectively.      For  a    compressor  blade    SWEEP   is  not   used. 

DIHED1,DIHED2 , DIHED  are  the  dihedral  angles  of  the  singular  line,  at 
the  hub,  at  the  cowling  (or  the  tip)  and  at  the  far 
field  respectively.  In  the  compressor  case  DIHED  is  not 
used . 
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Title  Card   5: 

Z,    XLE,    YLE,    CHORD,    THICK,    ALPHA,    FSEC 


Z 

XLE.YLE 

CHORD 

THICK 

ALPHA 


FSEC 


is  the    location  of   the   span   section, 
are  the   coordinates    of   the    leading   edge, 
is   the   chord   of   this   section,   used    to   scale   the   profile, 
multiplies      each  y   coordinate;    thus   the   thickness   of   the 
section   is   multiplied  by  THICK. 

is  the  angle  of  twist  through  which  the  section  is 
rotated.  Changing  the  value  of  ALPHA  by  the  same  amount 
at  each  span  station  introduces  a  stagger  angle, 
indicates  whether  section  coordinates  are  to  be  read  in 
from  data  cards.  FSEC  =  0  means  no  further  cards  are  to 
be  read  for  this  section.  FSEC  =  1  means  the 
coordinates  of  the  section  will  be  determined  by  the 
next   input    cards. 


Title  Card  6: 
YSYM,  NU,  NL 
YSYM 


NU,NL 


indicates   whether   the   profile   is   symmetric  or  not.      YSYM 

=   1    for  a    symmetric    profile  and  YSYM   =  0  otherwise.      For 

YSYM   =   1    only   the   coordinates    of   the   upper      surface      are 

read. 

are  the  number  of  points   on  the   upper  and   lower  surfaces 

respectively. 


Title  Card  7: 

TRAIL,   SLOPT,    XSING,   YSING 


TRAIL 


SLOPT 


XSING.YSINT, 
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is   the   included    trailing   edge  angle   in   degrees,      or      the 
angle      between      the      upper      and      lower     surfaces      if    the 
profile   is    open. 

is  the  slope  of  the  mean  camber  line  at  the  trailing 
edge.  SLOPT  is  used  not  only  for  the  profile  but  also 
to  determine  the  vortex  sheet  behind  the  profile.  If 
SLOPT  is  too  big  and  the  gap-to-chord  ratio  is  small  the 
program  may  abort. 

are  the  coordinates  of  the  singular  point  inside  the 
profile  about  which  the  profile  is  unwrapped.  This 
point  must  be  chosen  in  such  a  way  that  the  mapped 
coordinates  are  as  smooth  as  possible  without  a 
pronounced  bump  near  the  nose.  If  such  a  bump  occurs, 
we  first  move  the  bump  nearer  the  nose  by  positioning 
the  singular  point  closer  to  the  upper  or  lower  surface. 
Then  the  bump  is  smoothed  by  placing  the  singular  point 
at  the  right  distance  from  the  nose.  Moreover,  the 
singular  point  must  be  redefined  not  only  for  any  new 
profile  but  also  for  any  new  configuration,  since  the 
geometry   changes   with   the    parameter   DC. 


Title   Card   8: 

X,    Y 

X,Y 


are  the  coordinates  of  the  upper  surface.  There  are  NU 
cards  and  the  coordinates  are  read  from  the  leading  edge 
to    the    trailing   edge. 
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Title  Card   9: 

X,    Y 

X,Y  are      the    coordinates    of   the    lower  surface.      There  are  NL 

cards  and  the  coordinates  are  read  again  from  the 
leading  edge  to  the  trailing  edge.  The  first  data  card 
repeats  the  first  data  card  for  the  upper  surface  which 
defines   the   coordinates    of   the    leading  edge. 


These  cards  are  followed  by  (NC-1)  series  of  cards  beginning  with 
Title  Card  5. 

The  program  provides  both  graphical  and  printed  output.  The 
Calcomp    plots  were    shown   in   the    previous  section. 

In  the  printed  output  we  can  read  the  listing  of  all  the 
coordinates  given  in  the  input.  Then  the  mapped  surface  coordinates 
are  printed  for  the  hub  and  the  tip.  This  allows  us  to  ascertain 
whether  the  coordinates  of  the  singular  points  are  well  chosen.  After 
the  iteration  history,  with  the  maximum  correction  and  the  mean 
correction  to  the  velocity  potential  and  with  the  maximum  residual  and 
the  mean  residual  for  the  difference  equations  at  each  iteration,  the 
program  gives  us  a  printed  plot  of  the  coefficient  of  pressure  on  the 
surface.  This  is  followed  by  a  Mach  number  chart.  These  last  two  are 
given  for  each  span  station.  Finally  the  characteristics  of  the  blade 
are  printed. 

If     the     mesh  has   to  be  refined,  a   new  series   of   output   is   printed 
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for  the   new  mesh.      If      not,      the      program      terminates      or     restarts      by 
reading  a  new  data    deck. 
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VIII. 


LISTING    QF     THfc     CODE 


10 
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wk  IT, 430) 

v»RIT, 210) 

EAD,420)  TITLE 

wRIT,h60)  TITLE 

EAD,42G)  DESC 
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D  (7)  NX,NY,NZ,NM,K1,K2,J 

NX-H 

NY  +  2 

NZ  +  3 

30  K=1,MZ 

D  (7)  ( (G(I*J»K)« I=1»MX ), 

T  1NUE 

D  (7)  (E0(<) ,K*K1,K2) 

TINUb 


L  SU 

IL,  X 

C1»K 

(  IND 

(FCO 

1 

L  ES 

TINU 

=  J  I  T 

./SC 

TE  ( 

NT*C 

=  fit 

=  CQV 
A  =  BE 
NX  +  1 


KF  (ND,NE,NC,NX,NZ> IS 
s»Yi#NP»  ITE1,ITE2,IV> 
SYM  ) 

.EO.O)  GO  TO  190 
NT.Gh.l  .)  GO  TO  50 

TIM 

E 

AL 
IWRIT,A30) 

(NM)+JIT 

O(NM) 

TAO(NM) 


IT#U1»V1»W1 


J=1,MY) 


IP,XMAX,ZMAX,KSYM,PA,SY,SCAL,SCALZ,AX, 
,32,B3,Z»C1,C2,C3,C4,C5) 

CH0RD0»PA,SwEEPl,SWtEP2,SWtEP,DlHtDl,D 
»XZZ#YC»YZ»YZZ#Z,E1»E2,E3.EA*E5,IND»Z0 

YM,KTEl>KTE2>SCAL»A0>Z>ZS,XC>YC»SL0PTj 

S0>XP,YP»D1,D2»D3»X»Y,IND»XK*PA,>Z*YZ> 
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7C  C 


PO  C 


D 

W 

9C  C 


100 


110 


y  =  .\y 

7  =  Ul 
Y  =  \Y 
1=3 
?=NZ 
'lit 

0  7C 
kITE 
ONT  I 

»i  ie 

R  ITE 

1  6u 
RITE 
ONTI 
RITE 
P  iTt 
Rl  TF 
RITE 
0  90 
SI  IE 
'INT  1 
PITfc 
RITE 
^  I  T  E 
RITE 

0  10 
KITE 
RITE 
ONTI 
PITE 
PITE 
PI  IF 
RITE 
RITE 
RITE 
ALL 
PITE 
RITE 
PITE 
RI  TE 
X  =  NX 
L  =  0. 
3  il 
1=  I  T 
(  II) 

1  =  SC 
(  II  ) 
(LX) 
1-SC 
(LX  ) 
HO  PL 
ONTI 


+  2 
+  3 

♦  j 

+■^ 
( iw* 

1  =  2* 
(  Ik? 

NUc 
(IWF 
(  1  *P 
I  =  2» 
(  I  Wi= 

NOt 
(IwR 
(  i^P 
(  I  WR 
(  UP 
J  =  ^» 
(IWR 

N  UE 
(IwR 
(UP 
( IwR 

(IwR 

0  <=K 
ll.k 
(  UR 

NUE 
(  IWP 

(  UK 

(  IwR 

(  UP 
(IWR 
(  UP 

SEC  ON 
(  UP 
(  IWR 
(UP 
(  i  WK 

/2  +  1 


2  20  ) 

4  60)  (IV(I>K),K=K1,K2) 

4  30) 
2  30) 

440)  A0(  I)»50(l»Kl),Su(IUTE2) 


240) 
44  0) 

43C  ) 
250) 


2fcC) 

440) 

430) 

27C  ) 

2 

440) 

440) 


XMAX>  AX 


4  4  0)  8  0 ( J  ) 


SY»  AY 


Z(K)»XC(K)»YC(K)*XZ(K)»Y2(K)fXZZ('<)#yZZlK) 


280) 

440)  ZMAXt AZ 

430) 

2  90) 

3  0  0) 

47C)  NXjNY,NZ 
TIME  ) 
510)  TIKE 
310) 

4  4C)  FMACH»OM»AL 
320  ) 


0  K 
El  ( 
=  1. 
AL* 
=  xC 
=  1. 
AL* 

=  xC 
(<) 

NUt 


K2 


=  K1 

K  ) 

+  .5*SCAL*(A0(Il)*A0(Ii)-S0(Il»K)*S0(UfK)) 

AOdl  )*S0(  I1»K  ) 

<K)+AL0G(X(Il)**2*Yl**?)/XK/2. 

+  .5*SCAL+(AO(LX)*AO(LX)-S0(L<.l<)*iiO(LX,K)) 

AO(LX)*SC(LX,K ) 

(K)+ALQG(X(LX)**2+Yl**2)/XK/2. 
=  X (  I  l)-X(LX) 
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120 


K 

s 

D 

C 

s 
c 

A 

130  N 
P 
P 

P 

c 

L 

V 

w 

c 

R 
1 
R 
I 
1- 

W 

4 

C 
w 
L 
C 
I 
1 
J 

c 
c 
c 

1, 
I 

w 
W 
* 
w 
140    w 

I 
c 
c 

c 
1, 

V 

i 

w 

i> 

I 
I 

V 


ZDU*=KTb2 

=  0. 

G    120    K=K 

zo=.5*  (z( 

=  S  +  L  ZC* (C 
CNT iNUfc 
BLADE=S 
IT*Ml  +  i 
1=P10(NM) 
2  =  P  2  0  (  N  M ) 
3=P30(NM) 
ALL    BOUND 
2=U1+CA 
2=V1+SA 
RITE     (IWP 
DUNT(NIT) 
£S (NIT) =F 
F     (N1T.LT 
ATE=0. 
F     (NIT.GT 
COUNT (J  IT 
RITE     (ImR 
PITE     (IWP 
ALL     SECUN 
KITE     (IWK 
X=NX/2+l 
C    150     K=K 
F     (K.LT.K 
l»ITfcl(K) 
2=1TE2(K) 
AtL    VELD 
HORD(K )=X 
ALL    FORCF 
A3,C5(K), 
F     (KPLOT. 
RITE     (IWP 


-I 


TE1» KZDUM 

K+l )-Z( K)  ) 

H0R0(K  +  1  )+CHORD(K)  WC5(K  ) 


RITE 
RITE 
RITE 
RITE 


(Iwk 

(IWR 
(  IWP 
(  IWR 


150 


F  (KPLOT. 
ALL  SPEED 
0NT1NUE 
ALL  TOTFO 
PA,A8LA0E 
LD«0. 
F  (A8S(CD 
RITE  (IWP 
PITE  (UR 
F  (KSYM.E 
F  (KSYM.E 
=ATAN2( V2 


IT>490)  NIT,0G»IG,JG>KG»FR,IR,JR><R»GM,R*,F1,U2,V2,NS 

=NIT-1-JIT 

R 

.MIT  .AND.ABS(DG)  .GT.COV. AND. ABS(PG)  .LT.1C.  )  GO  TO  130 

.JIT+1)  RATE=(A8S(RES(NIT)/RES(JIT+1)))**(1./(CJUNT(NIT) 

+  1)  )  ) 

IT, 330) 

IT, 5  00)  RES(J1T+1),RES(NIT),C0UNT(NIT),RATE 

0  (TIME) 

IT, 310)  TIME 

1,MZ 

TF1.0R.K.GT.KTE2  )  GO  TO  150 


(k,K,SV,SM,CP»X,Y) 
(Il)-X(LX) 

(I1,I2,X,Y,CP,AL,CHQRD(K),XC(K),SCL(K),SCDK),SCM(«),A2 
PA) 

GT.1.AND.K.GT.KTE1)  GO  TO  140 
IT,  <»30) 
IT,34C) 

IT, 440)  FMACH,OM,AL 
IT, 350) 

IT,  4 40)  Z(K),SCL(K),SCD(K),SCM(K).CHORD(K) 
LE.l.)  CALL  CPLOT  ( 1 1, I  2, SM , AO, S 0 ( 1 ,K ) , FM ACH ) 

(K) 

R  (KTE1,KTE2,CH0RD,SCL,SCD,SCM,Z,XC»C5,CL,CD,CMP,CFP,CMY 

) 

J.GT.l.E-6)  VLD=CL/CD 

IT, 430) 

IT, 360) 

0.1)  U2«CA 

Q.l)  V2»SA 

,U2) 
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i 
L 
J 
D 
C 

IP 
160  C 
I 
N 
N 
N 
C 

1  A 

c 
II 

2N 

C 

IT 
2A 

I 
C 
N 
N 
G 
C 
w 
w 
D 
k 
C 
w 
E 
R 
& 
* 
u 
S 
C 

c 
s 


170 


1*0 


190 


200 


=  u;* 

1  =  1 . 
=  bGP 

PI  Tt 
RITE 

RITE 

sue 

RITE 

PLOT 

F     (L 

O  36 

ALL 

SC  AL 

ONTI 

F     (F 

X=  NX 

Y=NY 

Z  =  N2 

ALL 

Yi  AZ 

ALL 

HEC2 

E) 

ALL 

RAIL 

1,C1 

F     (  1 

ALL 

M  =  NM 

IT  =  C 

0    TC 

ONTI 

SI  TE 

RITE 

0    18 

RITE 

ONTI 

RITE 

NO     F 

EWIN 

0    TO 

RITE 

RITE 

TOP 

ONTI 

ALL 

TOP 


Aj-AL 
U2  +  V2 
/  F  M  A  C 

r/(Q/( 

( I  w  p 

(  I  UP 
(IrfG 

(  IwR 
<  iwR 
=  F*3L0 
PLOT. 
O./XH 
THREE 
) 

NUE 
HALF  ( 
+  NX 
+  NY 
+  NZ 
CCORD 
t  TY  iF 
S  INGL 
»DIHE 

SU^P 
,  xS,  Y 

,KS  Yvi 
N  D  .  t  C 
REF1N 
+  1 


+  V2  + 

Jl-0 
IT  ,4 
IT, 3 
I  T  #  <t 

IT, 3 
IT, 4 
T*(  1 
LT  .  1 
/PA  J 
D     (  I 


Wl* 
+  C. 
.2* 
40) 
70) 
40) 
60) 
AC) 
.-F 
)  & 
*ZC 
PLO 


Ml 
2 

C)  ) 
FMACH#QM»  AL#U#  V 

CL,CD, VLO 

CMP»CMR,CMY*ABLADE 
H  A  L  F  (  N  M  )  ) 

0  TO  160 

NF/CHukDO 

T,$V,S*,CP,X,Y,TITLf,DC,AL,ZONE,U,V,CHORGu>X:>CAL> 


>|M).tO.C.)  GO  TO  170 


(NX,NY,NZ,XTEu,ZTIP»XMAX,ZMAX,K$YM,PA,SY,SCAL,SCALZ,AX, 
,AJ,A1»A2»A3,8  0»B1»32»63>Z>C1>C2»C3,C4,C5) 

(NC,NZ,KTEl,KTE2,CH0RD0,PA,SwEEPi,SWEEP2,Skr.tP,DIhtrJl,r> 
0,ZS,XLE,YLE,XC,XZ,XZZ,YC,YZ,YZZ,Z,E1,F2,E3,64,E5,IND,ZG 


NE,NC,NX,NZ,ISYM,«TE1»KTE2»SCAL,A0,Z>ZS,XC>YC>SL0?T> 
fITEl»ITE2#IV*S0*XP,YP#Dl»02*D3#X,Y»lND»XK,PA,XZ»YZ» 


(ND, 
St  N? 
) 

.  C)  GC  TG  190 


\|X,NY»NZ>N1,K1,K2,NIT.U1,V1>W1 
90) 

((G(I,J,K),I  =  l,fMX),J=l,MY) 

( EG( K )>K«K1,K2 ) 


tO 
NUE 

(8) 

(6/3 
0  K=l 

(6  ) 
NUE 

(8) 
1LE  S 
D  8 

10 

(  I  *  R  I T  ,  't  3  0  ) 

(  I  WRIT, 400) 
0102 
NUE 

PLOT  (0.,0.,999) 
0101 


210  FORMAT  (20H0NYU  COMPRESSOR  CODE) 

220  FORMAT  ( 49H0 IN D  IC A T 1  ON  OF  LOCATION  OF  BLADE  AND  VORTEX  iHEET,27-i  1 
IN  COORDINATE  PLANE  Y  =  0 . / 27H0 ( <  I  V (  I  t K ) ,K=Ki , K2 ) , I = 2 , NX  )  ) 
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23C  FORMAT  (49H0CHGRDwISE  CELL  DISTRIBUTION  IN  TRANSFQRi 
1NU  MAPPtc  SURFACE  COORDINATES  AT  HU3  ANO  TJ.P/J5HO 
2\>  HUB  PRGFILE,15H    TIP  PROFILE  ) 

24C  FGn,MAT  (15HG   TE  LOCATION'  >i5H     POWER  LAW   ) 

(^mON  jfmal  cell  distribution  IN  TRANSFORMED 
) 

(1SHO       SCALE     FACTO-U15H  POWER    LAw        ) 

(  4dHOSjANwI  Sb     CELL     DISTRIBUTION     AND    SINGULAR 


2  5 C  FORMAT 

1    Y 
2  60  FORMAT 
270  FORMAT 
1   Z 
2,15H 
280  FOP-MAT 
F  0  P  I-  A  T 
FOKMAT 
FORMAT 
FORMAT 


>       »  1 5 


PLANE/15  H-J 


LINE/13HC 
XZ 
) 


290 
300 
310 
3  2  C 


,15H 

,15H   AN& 
I  ,  4H   J 
COR.tlOH 

SONIC  PTS) 

2,15H 

MACH  NO 


NZ       ) 
OF  A  T  T  A  C  <) 
»  4H    K   »l3H 
i  t  A  N  R  E  S  .  t  1  0 


wORK 


15H 


GM 


,15H        L  /  j 
YAw,  9X6HAEL ADE  ) 


QMEG 
) 
) 


,  i&H       X  SING    ,  15H       Y  SING    ,  15H 
yz      ,15H         XZZ      »13H         YZZ 
(15H0   TIP  L0CATI0N»15H     POWER  LAW   ) 
( x9H0ITERATIVE  SOLUTION) 
(13  HO         NX       ,  15H  NY 

(lpHO     MACH  NO    ,15H       OMEGA 
(  1CH0IT  ERAT  I0N,15H     CORRECTION  ,  4H 
1      PESiCUAL    »4H   I  >1H   J  >4H   K  ,10H  ^1 F  A  Ni 
2H  REL  FC1  1»10H    XSPEED  ,1GH   YSPEEG   ,10H 
330  F  0  k  M  A  T  (  1  5  H  0  MAX  R  E  S  I  D  A  L  1  *  1 5  H   MAX  RESIOAL 

15H   REDUCTN/CYCLE  ) 
340  FORMAT  (24HCSECTI0N  C HAR AC Tt R I  ST IC S / 1 5H0 
,15H   ANG  OF  ATTACK) 

(/13H  SPAN  STA1  I0N,12X2HCL,13X2HCD,13X2hCM,10X5HCHGRD) 
(22H0dLAUE  CH AR AC T ER I  ST  I C S/ 1 5HG     MACH  NO  1  ,15H 
15H   ANG  OF  ATTACK, 15H      MACH  NO  2  ,15H   UTV.  ANGLE 
(lb  HO        CL      ,  15H      CO 
(/ZX6HCM  PITCH,  6X7HC.1  RGLL,9X6HCM 
(IX, 1  4  H WRITE  ON  TAPE6) 
(24H0BAD  OATA, SPLINE  FAILURE) 
(8F10.6) 
( 10A8) 
(  1H1  ) 

(F12.5, 7F15.5) 
(1X*8E15.5) 
(1HC1JA9) 
(18.7115) 
(  1>  ,3ZIA) 

111D,E15.5,3I4,E15.5,3I4,2E10.3,3F10.5,I1C) 
(  c  E  1  5  .  4  ,  2  F  1  5  .  A  ) 

( 15H0C JMPUTING  TIME,F10.3, 10H    SECONDS) 

</5x3hFNX,llX3hFNY,ilX3HFNZ,10X5HFPLCT,9X5HXSCAL,9x5HP$CAL, 
19X5hGPlUX,9x5HGRIDY/lX,dE14.3) 
5  3C  FORMAT  (/4X7HFIT<NM),ex8HCGvO(NM),8X7HP10(NM),3X7HP2C(NM),dx7HP?G( 

1NM) #6X9HBETA0(NM) , 6X9HFHALF ( NM  )  ) 
540  FORMAT  (/5X5HFMACH,11X2H0M,14X2HAL, 12X3H  PA,12X4H  DC  /1X>5E15.5) 
END 


: 

IEGA 

350 

FOR  MAT 

3bC 

FORMAT 

I  A      , 

370 

FORMAT 

360 

FORMAT 

390 

FORMAT 

400 

FORMAT 

410 

FORMAT 

42  0 

FORMAT 

430 

FORMA  r 

440 

FORMAT 

45C 

FORMAT 

460 

F3RM/T 

470 

FORMAT 

48C 

FORMA! 

490 

FORMAT 

5  0C 

FOPMAT 

51  C 

FORMAT 

520 

FORMAT 
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1C 


2  0 


SUB 
11.5 
2PA, 
GEG 
DIM 
16)  , 
2  C  ( 
DIM 
IRE 
1WR 
RAD 
REA 
REA 
*RI 
KSY 
IF 
NC  = 
SWE 
SWE 
SWE 
DIH 
DIH 
DIH 

51  = 

52  = 
Tl  = 
T2  = 
ISY 
XTE 
CHQ 
ZCN 
K  =  l 
PEA 
PEA 
WPI 
ALP 
2T( 
ZS( 
IF 
IF 
ZON 
ZS( 
ZT( 
CDN 
PEA 
PEA 
WRI 
NU* 
NL  = 

IF 
N  =  N 
REA 


EGM     (ND,NC,NR.ZS,XS,YS»XLE»YLt»:>LUPT,TUAlL,*P, 
EP,DIHED1,DHED2,DIHED,kSYM,XTE0,CH0RDC,ZTIP,I 


ur  ine 

EP2, S  * 

D) 

TRIG  JEFIMTIQN  CF  6LAuE 

SIDN  X  S ( N D  > 1  )  *  Y  $  (  N  0  »  1 )  ,  Z  S (  1  )  >  X  L  t  (  1  ) ,  YLE(l), 

RAIL(l),  XP(1),  YP(1).  NP(1),  X(6),  Y(b),  Z(6)» 


SlUPT 
XC(c) 


YP,SwEEP 
SYMG,  Xl<  , 


(  i  ),  7T( 
»  YC(fc), 


SIGN  DESC (10),  A(o),  B(6),  D(6),  E(6) 
=  5 

=  D 

.29b  7fl 

I  R  E  A  J  ,  2  3  C  )  DESC 

IREAD,22C)  ZSYM,FNC,SwEEPl,SWEEP2,SWEEP,DIHEDl,DlHED 

(I  wR  IT, 270)  ZS YM,FNC,SWEEP1, SWEEP  2, SWEEP, DIHFD1,DI HE 

SYM 

C.LT.3. )  RETURN 

=  b*t  £P 1/ k AD 

=bwEEP2/PAP 

iWfcEP/RAC 

sQIHEDl/ RAD 

=D1HEU2/ PAD 

L1HED/RAD 

( S  W  E  £  P 1 ) 

(S*EEP2) 

(DIHED1) 

(D1HED2) 

1 


2,0IHEL 
D2  , DIHED 


D  (  IR 
D  (  IR 
TE  (  1 
HA=  AL 
K  )  =  ZT 
K  )  =  AL 
(FSEC 
(K  .GT 
E«(l. 
1  )  =  (  1 
1)  =  1. 
T1NUE 
D  (  IR 
D  (  IR 
TE  (  1 
FNU 
FNL 
(  YSYM 
U  +  NL- 
D  (  IR 


EA0,230)  DESC 

£AD, 22  0  )  ZT(K),XL,YL,CHORD,THICK»AL.FSEC 

hRlT,2  80)  ZT(K), XL, YL, CHORD, THICK, AL,FSEC 

/RAD 

(K)/ZQNE 

UG(ZT(K)+(1.-PA))*PA+(1.-PA)*ZT(<) 

•EQ.O.)  GO  TG  3  0 

.1)  GG  TO  20 

-PA)+PA*ZT  (1) 

,-PA)*ZT(l ) 


hA0,230)  DESC 
EAD,22C)  YSYM,FNU,FNL 
wPIT,290)  YSYM,FNU,FNL,ZONE 


.EQ. 1.  )  NL  =  NU 

1 

tA0,230)  DESC 
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3C 


40 


50 


£  C 
70 


PC 


90 


wRI 

kEA 

IiiR  I 

DG 

Rb'A 

UPI 

CON 

L=N 

IF 

PFA 

W&I 

DC 

REA 

WRI 

J  =  L 

XP( 

YP( 

CON 

GO 

J  =  L 

DO 

J  =  J 

XP  ( 

YP< 

CON 

W»I 

WftI 

WRI 

CON 

SCA 

C(K 

XLE 

YLE 

XX  = 

yy  = 

CA  = 

SA  = 

XC  ( 

YC( 

DO 

XS( 

YS( 

CON 

SLO 

TRA 

NP( 

CHO 

IF 

wRI 

WBI 

K=K 


L   ( 

TF 

u  ( 

Tt 

3C 

D  ( 

TF 

TIN 

L  +  l 

(YS 

U  ( 

Tc 

4C 

L  ( 

TF 

-i 

J)  = 

J)  = 

TH 

TO 


1 R  t  A  0 ,  2  2  C  )  TRL,SLT,XSING,YSING 

(1«RIT,3C0)  TRL,3LT,XSING,YSING 

1 R t A  0 1 Z 3  G )     DESl 

(IwRIT, 310) 

i=NL  »N 

U  E  A  0  *  2 .2  0  )     x  P  (  I  ) .  Y  P  (  I  ) 

(  I  w  R  I T  ,  2  6  0  )     X  P  (  I  )  .  Y  P  ( I  ) 

Ufc 

Y.^.GT  .0.  )     GO    TO    50 

I R E A 0 » 2 3  0  )     DESC 

(  IV/PI  T,  3  20) 

I»l»NL 

1READ»220)     V  A  L  »  C  U  M 

(Ikfll,2fc0)     VAL.DUM 

VAL 
DUh 
UE 
70 


I  =  iM  L  »  N 


60 
-1 
J  )  -XP  II 

J )=-YP( 

TINUt: 
T  E  IUF 
TE  iltit 
Tt  (IWR 
TINGE 
L[  =CHb« 
)=EXP(- 
(K )=XL/ 
(K)=YL/ 
XP (NL  )  + 
YP(NL  )  + 
CGSULP 
SIN  (ALP 
K)=SCAL 
K)=SCAL 
90  1  =  1, 
I*K  )  =  SC 
I,K)=5C 
TINUt 
PT(K ) =T 
IL  U  )=A 
K)  =N 
RD0=AMA 
( YSYrt.L 
T  T  (UP. 
Tt  (Iv-F 
+  1 


) 

I  ) 

IT, 200) 
IT. 250) 
IT, 240) 


ZS(K  ) 

TPL,SLT,  XiING» YSING 


0/(  X 

PA*Z 

ZONE 

ZONE 

(  XSI 

(  YSI 

HA) 

HA  ) 

£*(X 

E*(  Y 

N 

ALE* 

ALE* 


P(l  )-xP(  >t  L  >  )  /ZONE 
S(K)  ) 

+(XSING-XP(NL))*C(K)*THICK*SCALt 
+(YSING-YP(NL))*C(K)*THICK*SCALE 
NG-XP(NL ) )*C (K)*THICK 
NG-YP (NL  )  )*C  <K)*THICK 


X*C A*YY*rHICK*SA  ) 
Y*THICK*CA-XX*SA) 

(XP(I)*CA>THICK*YP(1)*SA) 
(THICK*YP(I)*CA-XP(I)*SA) 


AN( ATAN( TH ICK+SLT) -ALPHA )*C(K  ) 
TAN(THICK*TAN(TRL/RAD)*C(K)) 

XI (CHORDO, CHORD) 

E.O. . OR. ALPHA.  NE.O.)  ISYMO=0 

IT, 210)  ZS(K) 

IT, 250)  XL, YL, CHORD, THICK. AL 
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13C 


14C 


IF 

XK  = 
II 

DO 
CO 
x  (K 
Y(K 
Z  (K 
CON 
El  = 
E2  = 
Fl  = 
F?  = 
El  = 
E?  = 
CAL 
CAL 
CAL 
CAL 
DO 
XS( 
YS( 
CON 
CON 
DO 
X  (  K 
Y(K 
Z(K 
CON 
F  A  = 
EB  = 
EC  = 
ED  = 
El  = 
E2  = 
Fl  = 
F2  = 
CAL 
CAL 
CAL 
CAL 
DC 
X(K 
Y(K 
Z(K 
CON 
EA  = 
E6  = 
EC  = 
ED  = 
El  = 
E?  = 


(  «  .  L 
XK  *  ( 
(PA. 
120 
1  0  ) 
)  =  XS 
)  =  AT 

)=AL 

T  i  N  0 
1.  /  ( 
1./  ( 
(COS 
(CDS 
S1*E 
52  *t 
L  bP 
L  IN 
L  SP 
L  IN 
110 
1»K) 
I*K  ) 
TINO 
TINU 
13  j 
)  =  XC 
)=AT 
)  =  AL 
TINO 

cost 

SIN( 
COS  ( 
SIN( 

Sl/( 

S2*Z 

(EA* 

(EC* 

L  S  P 

L  IN 

L  SP 

L  IN 

140 

)=XL 

)  =  AT 

)  =AL 

TINO 

COS( 

SIN( 

CCS( 

SIN( 

Sl/( 

S2*Z 


{  i  .- 
N  e  .  1 
1  =  1, 
<  =  !> 
(  I,K 
AN2  ( 
C6(  Y 
E 

CuS( 
CLS  ( 
(Yd 
(  Y(  M 
1 

2*ZT 
LiF 

r  pl 

LIF 

TPL 

K»1j 

=  D(  K 

=  E(K 

E 

E 

K=l, 

(K) 

AN2( 

0G(  Y 

E 

Yd  ) 

Yd  ) 

y(n: 

Y(NC 

tA  +  T 

r  (nc 

Tl-E 

T2-E 

LIF 

TPL 

LIF 

TPL 

K  =  l» 

E  IK) 

AN2( 

0G(  Y 

b 

Yd  ) 

Yd) 

Y(NC 

Y(NC 

EA+E 

T(NC 


)  GO  TO  10 

pa) /CHCRDO+^A) 

.  0)  GC  70  150 

N 

NC 

) 

rS( I»K) »ZT(K) ) 

S  (  I*  K ) **2  +  2T (K ) **2  )  II. 

Y(1))+T1*SIN(Y(1))) 
Y(NC))  +  T2*SIiM(Y(NC))) 
>)*T1-SIN(Y(1)))+E1 
C ) )*T2-SIN( Y(NC ) ) )*E2 

(NC  ) 

(l»NC»Z>X,A,3,C»i#El»l>E2,C»G.tlND) 

(1»NC>ZS>D»1»NC»Z>X,A,?,C»0) 

d,NC»Z»Y.A.6»C»i*Fl,ltF2.C»0.*IND> 

(l>NCjZS#Etl»NC#Z#Y,A>B»C»0) 

NC 

) 

) 


NC 


Y  C  (  K  )  ,  Z  T  (  K  )  ) 

C  (<  )**2  +  ZT(K)  **2)  /2. 


) 

) 

)  ) 
)  ) 
1*5 
)/( 
b)/ 
0)/ 
(  1, 
(  1, 
(  1, 
(  1, 
NC 


B) 

EC+T2*EC) 
(  E  A  +  Tl  +  E  ti ) 
(EC+T2+E0) 

NC,Z,X,A,3>D,1,E1,1>E2,0,G.,IND) 
NC»ZS#XC»l»NC»Z>X»A»8tD»0) 
NC»Z>Y,A,3,D,l»Fl»l>F2,t>0.,lNL)) 
NC»Zj#YC*l»NC»Z#Y»A,B»D*0) 


YLE (K)»ZT(K)) 

LE( K)**2+ZT(K)**2)/2. 


) 

)  ) 

)  ) 

B*T1) 

)/{ EC+ED*T2) 
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15C 


160 
17  C 


IPC 


19C 


F 1  =  ( f. 
F2=  (F 

CALL 

CALL 

CALL 

CALL 

GO  TO 

CnNTI 

Z  0  =  Z  S 

co  u 

ZS(K) 

CGNT  I 
CONT  1 
CO  19 
DO  It 
XS(  I* 
Y  c-  (  I , 
CUNT  I 
XX5  =  F 
X  X  =  5  0 
X  X  r>  =  I 
X  T  E  C  = 
C0NT1 
ZTIP  = 
RE  Tlk 


A*  [ 

:*t 

SPL 
INT 
S^L 
INT 
1  I 
NUE 
(  I  ) 

J   K 

=  ZS 

r.uE 

NUt 

C   N 

0  1 
K  )  = 
K  )  = 
NUE 

xp( 

K.T( 

XX  + 

A'VA 
NUt 
Z3  ( 


i-cr )/ (EA+E3+T1) 

<:--_>)/  (  EC  +E  n*  T  ^  ) 

If  (lfNC»2#X*A#fl»C*l#61»l*t'2#0»0.»lN 

PL  (1»NC>ZS»XIE»1»NC»Z»X,A,8»D,0) 

IF  (l^C»Z»Y»A,-i,J»i.»Fl»i»FZ»0»U.»lN' 
^  L  (l»NC»ZS>YLE»l>NC»Z»r»A»i»L).u) 


=  i  »NC 
(  K  )  -  Z  0 


=  1  ♦  N  C 

■1»N 

Xi ( I »K )-XC ( K) 

YS(I*<  l-YC  U  ) 

XK*Xi(l>K))*CDS(X"<*YS(l#K))-l 
bX3(2.*XK*XS(i»f<))-i:.*XXS-l.) 
XXS)  /?. 

xl(XTEt»XXS) 

NC  ) 


[  ) 


200 
210 


220 
230 
240 
2^0 
260 
270 

260 


>F10.5-/15HO 

If  SING 


1H 

2PHA 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


) 

Z    »    t  F 1  o  .  5 
.13HTH1C* 


FORMAT     (16H0PR0FILE     AT    Z 
lLPt       ,  1  5  H  «    ilNG  » 1  5 ri 

FORMAT     (27HCSECTIQN    DfFINlTION     AT 
YLE  *lbH  ChORD 

) 

(8F10.6  ) 

(  1GA6  ) 

(  1H1) 

( F12 .4.  7F15.4  ) 

(8f 15.5) 

(/5X<tHZSYMjl2X3HFNC>10X6HSWtEPl>*X6HS 
HHEDl»9XeHDlHED2#lCX5HJlHED/lX»8ElD.t) 

FORMAT    (/5X5HZS(K  )  * 12X2HXL#13X2HYL» 11X5HCH0R 
12X<.HFStC/lX,7E15.5 ) 

(/6X<tHYSYM»llX3HFNU»12X3HFNL/lX,<tti:>. 

(/t>X3HTRL»12X3HSLT#HX5HXSIN3>10X5HYS 

(/5X5HXP(  I  )»1QX5HY?(I  )  ) 

(/6X3HVAL>12X3HDUM) 


ANGLE 

1 1  b  H  C 
NtSS    R 


XL 

r  i  o .  i  5 


U     >L 


.   1  3 

AL 


290 

FORMAT 

300 

FORMAT 

310 

FORMAT 

32C 

FORMAT 

END 

«tER2» 

0  ,  1 0  *  5 

j) 
ING/lXj 


9  X  5  H  S  *  E 
H-ThiC** 


EPjlGX'bniJ 
12X.2HAL  »  1 


*  1 1 5  .  S  ) 
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1C 


117, 

bt  r 

DIM 
1Z(1 

Q*  = 
DY  = 

K  Y  = 

Z0  = 

K?» 
C7  = 

K  1  = 
AEi  = 

ax  = 

AY=^ 

A/  = 

Hi  = 

*,WA 

IF 

373 

ZM,A 

IF 

SCA 

3  0  A 

T=2 

SY  = 

ax  = 

BX= 

V?  = 

hl» 

w?« 

bftX 
A^X 
CriX 
AP3 
MX  = 
LX  = 
LQ 
DD  = 
IF 
IF 
B»l 
IF 
A  =  C 
AS  = 
C  =  A 
DO  = 
01  = 
D2  = 
GO 
IF 
A  =  A 


t-  iMj  T  I  n  L     Z  J  J  R  0     (NXfNYpNZ*XTtO*Zl  I  ^*  X  I  AX  ,  Z  v.  A  X  *  *  S  Y  *  ,  y  A,  3  Y  ,  iC  AL  *  S  C  A 

A  X  ,  A  Y  ,  A  /.  *lY*F,AC»Al,Ak:.Aj»BC*ei/B2jft3»^»Cl*C?»C3fL^*Ct)) 

S     ;>H     3IKFTOEL     AM)     jPAN*13E     CDLSDIMTES 

h  N  S  1  J  M     4  C  (  1  )  ,     Aid).     A  >  ( i )  ,     A  3  (  1  )  »     B  0  (  1  ) »     r  1  (  1  )  ,     B  2  ( 1  )  ,     ti  3  {  i  )  * 

)  »     CHI)..     C  2  (  1  )  »     C  3  (  i  )  ,     l  *  (  1 )  *    C  !>  ( 1  ) 

2.  /NX 

1  ./Mr 

\Y+  1 

N  2  +  I 
1  .  /  i*  Z 

i  . 


X=  . 
(*3 

=  5  J 
x  =  . 

(  Kb 
L  =  X 

L  2  = 
.  IS 

S  U  K 
Art  i 
(  F  * 
ID* 
1  ./ 
Ul 
=-8 
=  1. 
=  (  1 

X=  A 

NX  + 
M  / 
3C 
(  I- 
(  1. 
(  1  . 


.  F  J .  1 .  )     XMA<*0.bZb 
(7  j  .  ) 


.  E  Q  .  0  ) 

0  /  (  .  5  J 
IP/ (  1. 
L 

1  )  /TY 
(Oi>?5 
/  3  -  1 1 . 
Y)  **2 
AL  I 

X/DZ  )* 
SO^iT  (? 
b  X  *  3  CM 
+  S  7  3  )  * 
+BdX* ( 


Z  M  A  X  =  1  .  -  i  .  /  N  Z 
GC1*XNAX*XMAX) 
C0C001*7MAX) 


0.*( T-0.U7)* Aob (T-0.C7)  ) 
)/?. +BX  »/ (F  +  l.  I 


*2 

.  +  ( 7.  +  S7  3  )  )/((l.+S73)*XMAX**3) 

T((7.  +  37j)/12.)*X,MAX**3 

X*AX*XMAX/i2. 
3.*C3X-<..*X^AX*Xi>UX)*XMAX*Xi'iAX/bOKT(C6X-X^AX*A,',AX) 


1 

1,  ^X 

)*  jx 

. 1 )  DO=-l.+l . /NX 

.NX)  0  D  =  1  .  - 1  .  /  M  a 


(ArjiiuD).Gr.xfAx)  c-u  r  o  i o 

ex- DO* DO 

SOhiT  (  A) 

Bx*Ai4btJX*(3.*CBX-4,*DD*DD)*DD*D9 

ABX*DD+dBX*AS*DD**3 

A3/C 

dexMUUM-6.*CttX+lSi.*DD*Di))-12.*0D**<t)*Du/U*C) 

TG  2u 

(  L  J  .  L  I  .  J  .  )  H  =  -  1  . 

d- ( (OC-d*XMAX)/(l.-XMAX))**2 


-  7, 


+  ; .  )  /  (  (  ■■     <■     )  ■   '      (   i  .  -  >  •  '  >  ) 


C  = /.  <  -  \  > 

I      --(/I    W     <,-,.)..(••-     O 

,         :        '■■.,*  +  .,        s  <    -I   (   i    ■-'-')  /, 

'   J  =  r  ;   -  /  (  (  .      <-,    I  <',<•-  ) 

[      =  -  l  .-  +  -:>):•(    i  „-••■:•*..>   ;  .   (   •  . 

;  ■■  (  i  )  -  . . 

r ;  ( I )  = .  -  i  .  ■ :  /  • .» 

(■    i  i  )  -    i  -    i 

,  =  •'.•-  ^ 

i     '  =  (  "   f  -  J  1       '  v 
:  -  .».  <  <  A  Y 

:-iuy4M-].)*u^-i) 
i;  1  =  a  *  C  /  (  U  -3  +  i  )  *  :  r  ) 

•      l  J  )  =  :  r  <■  J  'J  /  C 

H  2  i  0  )  =  L*  i  *  C 1  *  v  ? 

■  ;  (  .  )  =  -  4  Y  *  )  i  v  j  Y  *  (  3  .  *  i  ■  +  J  )  /  I  I  a  i  +  0  )  *  A  ) 

-  ->  Z  =  -  "^  *  i  C  -!  T  (  •>  .  *  (  7  .  +  i  7  '<)  )  /  (  I  1  .  ♦  b  7  3  )  *  Z  >  A  X  *  *  3  ) 

t  ;•'  L  =  1  .  -  i  i-  /  *     „ -  T  (  (  7  .  +  >  7  3  )  /  i  2  .  )  *  Z1"!  A  X  *  *  j 

=  (  is- .  +w  i  )  *  7  h  A  >  *  7  i"i  AX/  i,>., 
{. .-  Li  i  =  A  -i  L  +  J  o  £  *  I  2  .  *  C  p  7  -  *,  .  *  >  •/■  U  a  *  Z  (■  l>  *  •)  *  Z  •'  A  x  *  i  y- 1.  > 

;  -1    v  i     k  =  c  >  \  > 
:  ,:•=  (  k-m  )  1  J/-/, 
ic     (k  .t-..vZ)       .L^I")->.  *0Z 


I  11. 


-  4.  ?    -5  ' •  *  ^ 


■    ) 


it-     ( /  6  o  ( [.  i ) )  .  C  T  .  i  V-  a  X  )     u    T ,  J    jo 

A  =  Cti?-L-L)*l   i< 

&S=SC^T(a) 

C  =  A f .7  *  As  +  r  j  i:  * '(  3  .  *'C  -J  ?-"*  •.  *  L-  0* L D  )  *  J i>+ 1   ) 

P.0=  At  Z*  JU  +  );•  7  *  AS*Q-'D**  3 

C 1  =  A  S  /  C 

[)?=VjF  7*  (  C  J  Z*  (  -h  .  *Cu/+  1  -.  *  JU*L-D  )-!?.  *  J  j**^  )+  --  /  (  a  *:  ) 

GU     1-     oC 
:  (      IF     ( [ U . L 1 . 0  .  )     b  =  - 1  . 

A  =  i  .  -  (  I  u  b-  o  *  7  '.  A  *  )/(l.-Z^AX))**2 

C  =  A  *  *  A  Z 

D=( A/+AZ-1 . )*( 1 .-A) 

Lu=  B *  Z  •-!  A  \+  A 3  !i  Z  *  (  IJ D-B  *  Z  \*.  A  X  )  /  C 

Cl«A*C/({l.  +  L-)*Ae3Z) 

D?  =  -(AZ  +  A7)*(OL3-6*ZivIAx)*(3.+D)/{(l.+D)*A*{i.-7^AX)^*£;) 
^0     Z(K)=jCALZ*L)0 

C4J(K)  =  Ji*jl*k<' 

C-5(k)=c<P(-?;,A*7(K)) 
70    C^.(^)=CD('<)'*C?(X) 

Z(K2+1)=Z(K2)+Z(K?)-Z(K2-1) 

h  c  r  u  t-  \ 
i  no 
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1C 


20 


30 


1HK) 
2 1  NO 

G  E  N 
D  i  M 
1  )  . 
D  1 
EM 

L£>   ( 

COW 
Kl  = 
K2  = 
KTE 
DQ 
IF 
IF 
CON 

a  =  c 
si  = 

S2  = 

Ti  = 
T2  = 
1  F 
F  1  = 
F2  = 
Tl  = 
T2  = 
SI* 
S?  = 
CON 
CAL 
C  AL 
CAL 
C  AL 
CAL 
CAL 
CAL 
CAL 
b>  =  B 
Sl  = 
S'* 
T  =  3 
Tl  = 

r?  = 
>c  ( 

YC  ( 
N=K 


c  j  vJ  T  I  U 
1  #  MHtU 
, ZJNb  ) 
t  Y  A  T  t  i 
E  N  b  j  J  N 
Y  /  /.  (  I)  , 

1  C      K  =  1  , 
K  )  =  [j  . 
K)  =  j. 

7  I  MLb 

NZ  +  ^ 

1*<U 

2  C  K  =  *  1 
(  2  (  K  >  .  L 

(  Z  (  K  )  .  L 
T  1 N  U  E 
H&kDG/Z 
TAN(Swfc 

T  A  N  (  S  W  h 
T  A  N  (  C  I  H 
TAN ( Qlh 
(PA  .At. 

i .  /  ( c  a  s 

1./ (COS 
(CGi<  Yl 
(CCJ(YL 
Sl  +  Fl 
S2*F2*E 
T1NUE 


5  J  N  G 

ZS  (  1 

/  (  1 
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GL     ( NO  r  Z »k  rhA,K  TE2>CHGRD0* ^ A ,  S  *FF b j  ,  s,  ,  £3  2 
HE&jZS>XLOYLE»XOXZ>XZZ*YC>YZ>Y2Z>Z#F-J>E2> 

L  L  A  &     LINK     F  LI  *     TRANSFORATION 

)  >      X  L  r  (  1  )  .     Y  L  L  (  1  )  f      *  C  (  1  )  .     X  7  (  1  )  *     X  Z  Z  (  1  )  t      <> 

)  »     c  1  (  I  )  ,     t ;  2  (  i  )  ,     E  3  (  1  )  ,     E  h  (  1 )  i     b  5  (  1  ) 


»    .">  A  C   '     -    «    I.    J 

:  3  »  fc  *«  i  E  b  > 


(  i  )  .    y  ' 


T.  Z^ 
E  .  ZS 

DNE 

EP1) 

EP2  ) 

EDiJ 

ED  2) 

1.0) 

(YLE 

(YLE 

Ell) 

E  (  NC 


(  1)   )      KTE  i  =  K+  1 

(NO)    K  T  b  2  *  X. 


GG    TO    30 
(1))*T1*SIN(YLE(1)M 
(NC  )  )*T?*SliM(  YLE  (NC  )  )  ) 
)*T1-SIMYLE (i) ) )+Fi 
)  )*T2-SIN(YLE(NC )) )*F2 


XP  (Z  MNC  )  ) 


5  P  L  I F 

INTPL 

I K  r  pl 

INTPL 
SPLIF 

INTPL 
1 N  T  P  L 
INTPL 

*  T  A  N  (  S  U 
B*S1 
B*S2 

*  T  A  N  I  D I 

B*T2 
2  )=3.*( 
Z  )  =  3  .  *  ( 
T  E  2  + 1 


(  1* 

(  KT 

(K  r 
(K  r 

(  i» 

(KT 

(  KT 
(KT 

dbp ) 


NC*ZS»XLE»cl»E2*E3»l#Sl*l»S2»C#G.*lND) 

Fl,<Tt2»Z»XC.A,NC.ZS»XLF»tl,£2>t3.0 

E1j<TE2jZ»XZ*4.»NC»ZS*E1»E2»E3*E4»0) 

El#KTE2#Z»xZZ#l»NC#ZS»E2»E3jE4#t5#0) 

NO  Z  S»  YL  E >  El  >  E2>  E3 ,  1,  T  1»  1*  T  2,  G »  C  .»  INL ) 

El>KTE2»Z»YOi>NOZS»YLt»ti»E2*E3.0) 

El*KTE2»Z«YZfl>NC»ZSjEl»E2fE3»E*»C) 

El#KTE2»Z»YZZ»l»NC#ZStE2»E3*E^»E5#0) 


H&  0  ) 


XC  (3 
YC  (  3 


)-XC  (<»)  )  +  XC(  5  ) 
)-YC(4)  )+YC( 5) 
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40 

5C 


IF     (r-  .GT.K2)     GD     TO     50 

OG    4  C    k  =  N ,  <\  ? 

ZZ=(Z(K)-ZUTE2) )/3 

A=EXP(-ZZ) 

XC (K  )  =  XC(KTE2)+S*ZZ+{S2-S)*ll.-A) 

YC(K)  =  YC  UTh2)  +  T*ZZ+(T2-T)*ll.-A) 

XZ(K)=(S+(  j2-S  ) *A)  /3 

YZC<)  =  (  I  +(T?-T  )*A)  /fl 

XZZ(K)=-(S  <:-S  )*A/  (6*6  ) 

YZZ(K) =-(T2-T) *A/( b+6) 

CON! 1NUE 

CONTINUE 

YC( Kl-1 )=YC (Kid) 

YC ( K2d ) =YC ( K2-1 ) 

^  F  T  U»*  S 

END 
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2  0 


SUB 
1SL0 
2xZ, 
INT 
INT 
DIM 
1  XC 
?f    0 
31(1 
PI- 
SI- 
T  =  l 
D<  = 
LX- 
MX  = 

rz  = 
I  vo 

I  VI 

CO 

ITE 

ITE 

DO 

IV( 

S0( 

CON 

K=K 

K2- 

K?  = 


ROUTINE  S  U »  F  (NL),NE»NC»NA#NZ»ISYM»KTtl,«TE2»SCAL>AC>Z>Z:>>XC*YC» 

PT»TRAlL»XS»YS*NP#ITtl»ITE2#IV»S0»XP*YP>Dl>D2*D3*X,Y#lND»XK,PA, 

Y  Z  ,  A  1  ,  C  1  ,  <  S  Y  M  ) 

EkPOLATES  CAPPED  BLADE  SURFACE  AT  "Pbri     POINTS 

EPPQLATION  IS  LINEAR  FOR  CYLINDRICAL  COORDINATES 

ENilUN  39(nE»1)»  XS(ND.1)»  YS(NDd).  ZS(i)»  SL0FU1),  T«Ail_(l)> 

(l)i  Y  C  ( 1  )  >     A  0  (  i  )  ,  Z  (  i  )  ,  X  (  1  )  ,  Yd),  X  H  (1  >  ,  Y  P  (  1  )  ,  D 1  (  1 )  »     D  2  ( 1 ) 

3(1),  I  \l  (  N  E  i  1  )  t     N  P  ( 1 )  t     I  T  E  1  (  1  )  »  I  T  E  2  <  1  )  ,  X  Z  ( 1  )  ,     1 1  ( 1 )  t     Al(l),  C 

) 

3.14lt>92b5 

.i*SCAL 

./SI 

2.  /NX 

N  X  /  2  + 1 

NX  +  1 

NZ  +  3 

-1-ISYM-I  SYM-ISY1 

--1-ISYM 

1 C  K  =  1  f  1 Z 

1 (K )=MX 

Z(K)=MX 

10  1=1,  MX 

1*<  )=-2 

I,K)=0. 

TINUE 

TE1 

1 

K2+1 
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3  0 


VD 


bO 


7C 


80 


f  1 
R2 
U 
P2 
R2 
Rl 
C  = 
0  = 

ex 
o 
c- 
cc 

Dl 
IF 
If 

cc: 

IT 
IT 

CC 
KK 

Pa 

h  = 

XA 
YA 
KB 

YB 
Cx 
C* 
Q« 

CC 
X  ( 

en 

AN 

v» 

DO 
x  I 
YI 
R> 
IF 
AN 
U" 
V« 
R« 
XP 
YP 
GO 
AN 

U" 

v» 

xp 
YP 


«K2-1 
=  I . 

(ZS(K 
=  (EXP( 
=  PA*Ri! 
■1.-R2 
Pl*xS( 
R1*YS( 
=EXP( X 
i>  CJ(-  T  { t 
(C+Cx) 
=  S  C  R  T  ( 
5C  1  = 
(  (  AC( 
(  (  A0( 
NTIMUc 
El(K)i 
E2(«)= 
■AO( 12 
»K1 
Rl 

NP(KrO 

*EXP( X 

»EXP(  X 

=  ExP<  x 

»EXP(  x 

=SQRT( 

=  (XA+C 

SQPKC 

7C     1  = 

)  =  C*A 

TINUE 

L  =  PI  + 


2  )-Z  (<  )  )     20»4C#30 

Z(K))-EXP(2S(K1)J)/(EXP(ZS(K2))-EXP(Z5(K1))) 

+(l.-PA)*{Z(K)-ZS(Ki))/(ZS(K2)-ZS(<l)) 

1»K1 )+P2*XS<l»K2 ) 

1><]  >+P?*YS(l>«2) 

K*C )*CCS( XK*D)-1. 

XP (2. *  XK*C  )-2.*CX-l.  ) 

II  . 

(C  +  C  )/SCAL) 

2*'^X 

I  )*-.  5*DX  )  .LT.-CC  )     11*1+1 

I)-.5*DX)  .LT.CC)     12= I 


II 
12 
)/CC 


K*XS(1»<K))*C0S(XK*YS(1»K<))-1 

K*XS(l»KK))*SlN(XK*YS(l»KK)) 

K*XS(N>KK) )*CCS(XK*YS(N*KK) )-l 

K*XS(N#KK)  )*SIN(XK*YS<N#KK) ) 

XA*XA+ YA*YA  ) 

X)/2. 

X/C)  /CC 

1»MX 

0(  I) 


PI 


I  =  1  *  N 
P(XK*XS(I»KK))*CUS(XK*YS(I»KK))-1. 

P(XK*XS(I»KK))*SIN(XK*YS(I»<K)) 
T  (  XI**2  + Y  1**2  ) 

•LQ.O.)     GG    TO    80 
ANGL+ATAN2<(U*YI-V*XI)#(U*XI+V*YI)) 


T(  (R+R  WSCAL  ) 
=R*CDS( .5*ANGL) 
»P*SIN( ,5+ANGL) 

90 
PI 


I  )  *  0  . 
I)=0. 
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9  0  CON 
ANG 
ANG 
ANG 
ANG 
ANG 
ANG 
ANG 
Tl  = 
T2  = 
CAL 
CAL 

DEF 

X  A  = 

XI  = 

A=S 

e=i 

ANG 
0  =  1 

v  =  o 

M»I 

x  xs 

YV5 

YYS 

DO 

t=l 

100  CON 
E  =  l 
XX  1 
D=B 
YY1 
XX  = 
Y  V  = 
R»S 
ANG 
U=X 
V=Y 
P  =  S 
XP( 
IF 
Y(  I 

11C  CON 
X3  = 
A  =  S 
B  =  l 
ANG 
U  =  l 

v  =  o 

M  =  I 
XXS 
YYS 


TINUE 
L=  ATAN 
Ll= ATA 
L  2  =  A  T  A 
Li  =  ANG 
L5r  =  ANG 
Ll= ANG 
L2=ANG 
T  A  N  (  A  N 
T  A  N  (  A  N 
L  SPLI 
L  I  -»T  P 


(  SLO* 
N(  YA/ 
N  ('  Y  8  / 
L-.  b* 
L-.  5* 
ll  +  XK 
L  ?+XK 
GLl  ) 
GL2) 
F  (1, 
L  (II 


T(KK  )  ) 

X  A  ) 

XB) 

(ANGL1-TP  AIKKK)  ) 

( ANGL2  +  TkAIL ( KK )  ) 

*f3 (1»KK) 

*YS (Nf KK  ) 


N>XP*Y?Ol*d2»D3»l*Tl>l*T2#0»0.#IND) 
,I2>X,Y,1,N»XP,YP#D1»D2»D3»0) 


IMTIIN    OF     THE     VORTFx     SHEET 

XS(  1,KK  ) 

.25*XA 

LGPT(KK)*(XA-X1) 

./ ( XA-X1) 

L=PI+PI 


-1 
l.+S 

SCAL 
ATAN 
10  I 
/I  .0 
INUE 
05*E 
ALOG 
(XXI 
YYS  + 
XP  (  X 
XP(  X 
PT(X 
=  ANG 


CAL*(X(Il)**2-Y(Il)**2)/2./G/0 

*X(I1)*Y(11) 

2( YYS»  XXS) /XK 

■1»M 

5 


(l.  +  .5*E*SCAL*X(I)**2/0/ti)/XK 

-X  1) 

A* (0-1.  )/D 

K*XX1)*CQS  (XK*YY1)-1. 

K*<X1)*SIMXK*YY1) 

X**2+YY**2 ) 

L+ATAN2((U*YY-V*XX),(U*XX+V*YY)) 


KT  (  (R+*)/5CAL  )*U 

)=P*CDS( .5*ANGL) 

XP(1  ).oT. x(  I  )  )  GO 

*R*SIN(  .5*ANGL  ) 

INUE 

i ( N, KK) 

OPT( KK) *( XB-Xl  ) 

/( XB-Xl ) 

=  0. 


TO  100 


2  +  1 

=  1 .+SCAL*( < ( 12) **2-Y(  12 )**2  )/2. /Q/a 

»SCAL*X( I2)*Y(I2) 
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120 


13C 
14  0 

15C 

160 


170 


ieo 


190 


200 


YYS  = 
DO  1 
E  =  l  . 
CONT 
t=l  . 
XX1  = 
C=B* 
YY1  = 
XX=  E 
YY=E 
R-SQ 
ANGL 

u  =  xx 

V  =  YY 
P  =  SO 
XP(  I 
IF  ( 
Y(  I  ) 
CONT 
C=P* 
DG  1 
S0(  I 
CQNl 
IF  ( 
KK  =  K 
P=R2 
GO  T 
CO  1 
IV(  I 
CONT 
S0(1 
SO(M 
M=I1 
00  1 

ZZ  =  Z 

IF  ( 
CONT 
M  =  I2 
DO  1 

zz  =  z 

IF  ( 
CONT 
K2  =  K 
K=K  + 
IF  ( 
Kl=3 
K2  =  N 
DO  2 

zz  =  z 

IF  ( 
CONT 
K=K  + 


A  T  A.nI 
20  1 
/1.0 
1NUE 

05*e 
AlOG 
(  XXI 
YYS  + 
XP  (  X 
X  P(  X 

RT(  X 

=  ANG 


2 ( rYb, XXS  )  /XK 


(l.+.5*F*SCAL*X(I)**2/C/Q)/XK 

-XI) 

A*( D-l  .  )  /D 

K*ul)  *CCS  (XK*YY1 )-l. 

k*xxi ) *^  in( xK^rr i) 

x**2+  YY**2  ) 

L  +  ATAN2({U*Yf-V*xx),(U*XX+\,*YY>) 


KT( (k  +  X  >/5CAL  )*0 

)=h-*CC:>(  .5*ANGl) 

XP(  I  ).LT.X  (  I)  )  GO  TJ  120 

=  R  *  S  I N  (  .  5  *  A  N  G  L  ) 

INUE 

CC 

4C  I=1»MX 

>K) =S0( 1 »K )+Q* Y(  I  ) 

INUE 

KK  .EL.K2 )  GC  TO  150 

2 

0  60 
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INUt 

,K)=SQ(2,<)+S0(2,K)-S0(3,K) 
X*K)=SO(NX,K )+SO(NX*K)-SO(NX-l,K) 
-1 

70  1=2.^1 

(K  ) 

ZZ.GE.Z(KTE1  )  )  IV(  I,K)  =  1V0 

INUE 

+  1 

60  1=M#NX 

(K) 

ZZ.GE.Z(KTE1 ) )  IV ( I>K)=I VO 

IMUE 

2-1 

1 

K.Lt.KTE2)  GO  TO  20 

Z  +  2 

00  I=2»NX 

U) 

ZZ.LE.ZS(NC).ANO.ZZ.GE.Z(KTED)  IVC  I#K  J-  IVO 

INJE 
1 


7  ■■'  - 


21C 


220 


23  0 


IF 
N=K 
1  =  1 

01 

IF 
i  F 
IF 
C  0  N 
b  =  S 
IF 
CON 

n 

bO( 
S0( 
S0( 

CON 
1TE 
ITF 
ITfc 
ITF 
F  FT 
ENC 


(«  . 
U  I 

l  fi 

?c'J 

i  - 
(  ]  v 

(  IV 

(  1  V 

TI  M 
0(  L 
(?. 

1  IN 

2  30 
I  f  M 
It  * 
1»< 
T  IN 
1(  f< 
IK 
2(M 
?(< 
IF  N 


L  E  .  K  2  )     G  Li     T  U     iVu 


UTE 

1=2 

(  I  »K 
(  1  <-l 
(  i  +1 

X,K  ) 
/  K  r  * 
IE 

1  =  1 
Z)=3 

Z)  =K 
1-1) 

LIE 

Z)  =  K 
1-1) 
7  )=N 
1-1) 


1  ) 
If  K2 

,  N  > 

) .  c-  r . 
,K  +  i ) 

»  K  -  I  ) 


)     G  !     I  U    210 

G  T  .  0  .  I  F  .  I  >/  (  1  -  1  ,  K  +  1  )  .  G  T  .  0  )      1  V  (  I  j  K  )  =  I  V  1 

G7 . 0 . C* .  I V ( I- 1 >  *- 1 )  . G  T . G )     I  V  1 1  #  <  )  *  I  V  1 


5/(T-S*S).LT.l.E-u5)     W  (  L  x  ,  K  )  =  J 

,  MX 

,*(SG(I»M7-l)-S0(l»MZ-2))+SO(I»NZ) 
>YM*S0(I#MZ)+(1-kSYM)*S0{I>MZ-2) 
=  SC(  I>Kl  +  l) 

SYM*ITEl(MZ)  +  (l-Ki)YM)*ITEl(  M7-2) 
■ITEKK1  +  1) 
K  +  2-1TEKMZ) 
=ITE?(K1+1) 
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SUB 

INI 

COM 

I>  AG 

2,  CI 

3  Z  ( 1 

AHA. 

T  =  2 

KY  = 

^z  = 

on 

DQ 
DO 
G(  I 

COlNi 

KZ  = 

DO 

DO 

IF 

DSI 

DSK 


U  7  1 N  E 
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2  9  )  ,  A 
7)  >G2 
»  Y  Z  Z  ( 
>  SA,F 
SC  AL 
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1  =  1. 
J=1j 
K=l, 
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Not 


F^T  iv 
ri".ATE    OF    KEDUCEO    POTENTIAL 

2  9 »  2  t  *  1  7  )  »  S  C  (  1  c  *  *  i  7  )  »  E  0  (  1  7  ) »  I  V  (  I  2  -)  .  1 7  )  ,  1  T  E  1  (  i  7  )  ,  17  F  2  (  1  7  ) 
1  (  1  2  <? )  »  A  2  (  1  2  i  )  ,  A  3  ( 1  2  9  )  »  -3  I  (  2  6  )  »  M  (  2  6  )  ,?.  2  (  1  c  )  ,  i  3  (  2  fc  )  » * :  (  17  ) 
(17),C3(17),CM17)>C5(17)»XC(17),XZ(17),XZZ(17),YC(17),t 
17),KSYM,NX*NY*NZ><7E1»KTE2,ISYM.SCAL»SCALZ»XK,0MEGA,ALP 
MACH»T»PA,ABLAOE 
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It 

1  7 


1=2, NX 

K  =  K  Z»MZ 
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0(I+1»K )-SC( 1-1, K) 
J  (I>K+1  )-S0(  ItK-1 ) 
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I  X  = 

57  = 

*  =  A 

Y  =  5 
XX  = 

yx  = 

FH  = 
H=l 
A2  = 
8Z  = 
A  =  C 
Yl  = 

P  =  A 

Q  =  A 
F7Z 
FYZ 
F*7 
FYY 
FX  Y 
F  xx 
A*  = 
AY  = 
A7  = 
BY  = 

V  =  i 
U  =  C 
w  =  C 
G(  I 
CON 
CON 
Kl  = 
K?  = 
DG 
EC  ( 
CON 
PPT 
END 


Al  (  ]  ) 
CI  U  ) 
C(  I  ) 
0  <  I  .  s 

c  .  I  xK 
2  .  /  X  r. 

X  X  *  X  X 

. /Fh 

-nn 

-X  X*  Y 

5(K) 

YC(K) 

*C  DS  I 
♦  SIM 
=  A*  A 
=  H*FZ 
=  H*FZ 

r  (h*H* 

=  H*H* 
F>  Y-S 
F  YY-  S 
FYZ-i 
AY-S.X 
A*Xx* 
A*XX  + 
A*XA  ( 
»KY  +  1 
T  I  N  L  t 
T  1 N  l  E 
KTE1 
KTt2 
40  K  = 
K)=0. 
T  1 N  U  t 
URN 


*  D  S  I 

*  l)  S  < 


*X*<T*-X**?  +  Y**2)/{  (T+X**2-Y+*2  )**2+(2.*x*Y)**6) 
*  Y  *  ( T->**2-Y**2)/< (T+X**2-Y**2)**2+(2.+X*Y)**2) 
+  Y  X  *  Y  X 

ZU)-YX*YZ(K) 
A ( \ ) +YX*XZ  (K  ) 

+  ArAIM?(2.*X*Y>T+<  +  x-Y*Y)/XK 
PA  *Y1  ) 

PA^Yi ) 


Z*3Z 

i*AZ 

( YX  *Yx  +FZ 

(-XX*YX+F 

(  X  X*XX  +  FZ 
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X*F  X  Y-!>7* 
X*Fx  Z-S  Z* 
*AX-SZ* AZ 
F  /F77-C  A* 
SA*YX*P/F 
K ) +S  A* ( C  + 
»K )  =  G(  I»K 
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Z*( 6Z*riZ+  Xx*XX  )  ) 
ZZ*(AZ*8Z+XX*YX)) 

Z*( AZ*AZ+YX*YX ) ) 

FXZ 

FYZ 

FZZ 

YX-UMEGA*XX/C4(K) 

ZZ-OMEGA* YX/CMK ) 
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Y-l»K)+(AX*U+AY*V+AZ*W)/(BY*ei(KY)) 


SUBROUTINE     BOUND 

OFFINES  THE  BOUNDARY  VALUES  OF  THE  VELOCITY  POTENTIAL  G 
COMMON  G<lZ9»26fl7)*S0U29*17)*E0(17)*IV(129»17)»ITEK17)#IT£2<17) 
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2C 


COM 
1KG 

COM 
19)  , 
2,Q1 
LX  = 
MX  = 
KY  = 
PY  = 
MZ» 
DX» 
Tl- 
AAO 
ZSY 
Gl  = 
02* 
RM  = 
GM  = 
NM« 
FR* 
IR» 
JR«= 
KR  = 
DG  = 
IG* 
JG  = 
KG  = 
NS* 
II- 
12  = 
DO 
DO 
Gl( 
G2( 
CON 
JT- 
CON 
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KM  = 
KN  = 
I*L 
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DSK 
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DO 
G(l 
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IF 
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0 

0 

0 
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NX 
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TINUE 

1 

TINUE 
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t  =  r-*(T-rh'*>.J)**1i/(l.-K  +  {?./XK*YP)**£) 
,  /=-YPxT  /  ( •<  .*  )*p./X*/n-YP**2) 
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=  (H*(  .; 

k=  1+  I  ? 
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(  j  >  J  >  1 

r>- t  l  wu 

=  K¥u 
(  1  f  J  .  ". 

1  (  J  ,  J  ) 


1  +  1  ,  X-  M  )  -  S  C  (  1  -  i  >  K  K  ) 

i  »y  ;  )-i  r (  l  »kl  ) 
i  ( i  )  *  i  j  i 

L)=j(I»1»KM)4ZSYM*(G11»2*KL)-G«I*2»KN)J 


Y  =  H 


I  .  K 
(  T  + 
(  T- 
+  Y  x 
M  )- 
V!  )  4 
■    )* 

*  41 

*  5  ( 
) 

) 

r  )  - 

t  K  -i 

L  i  ( 

.  J  > 

.  J  t 

(  I  , 


r  ) 

)(**2'4Y**?)/((l4X**2-Y**2)+*i:+(2.*X*Y),)'*2) 

X*^-V**^l/(  (T4X<«*2-Y**?)**i:4(  2,*X*Y)**2) 

*  YX  ) 

Y  X  *  Y  Z  (  <  M  ) 

Jl 

(  I  ) 

(  I  )  ) -JT*-  Z (  1  )  ) *Hi  (  J  )* JT 


G  (  1  iv  »  J  »  k  N  ) 

)-G(I»J-li*M) 

KM)*G{I#J»KN)+A*{&{IF»J»KL)*0GI)-6*DLJ)/(C1(KM)4A) 

H  L) 

kn)+3.*(o( [»J»KL)-G(I»Ji<M) ) 

J  J4ZSYK* (G( 1 t J  » 1  )-G2(  I. J  )  ) 


l):(CH'<vi)*iJ(I>J>K\|)4AMG(  IP>JfKL)4UGl  )-6*DGJ)/(;il^)U) 
=  G  (  I  »  J  »  K  L  ) 
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7C 


IF 

CON 
IF 

KK  = 

L  =  3 
I  =  t> 
J=N 

K  =  L 
X  =  A 

Yl* 

DEN 
XX* 
yx  = 
H=l 
P  =  C 
G  =  C 
AZ  = 
BZ  = 
SX( 
SZ( 
GX  = 
GY  = 
GZ  = 
U  =  G 

V  =  G 
w  =  G 
00  = 
VO  = 
WO  = 
ul« 
vi  = 
wl  = 
OD 
IS* 
IT  = 
DO 

X  =  A 

Y  =  S 
DEN 
XX  = 
YX  = 
Yl* 
G?( 
Gl( 
G(I 
G(l 
SX( 
SZ( 
DSI 
DSK 
DSI 


(  1  .1  I  .L X  )  i 

1  I  Nut 

(  J  1  .  tu.-i  ) 

? 


TO  4  0 


GO  TO  liC 


&(  1 
G(  I 
P  A* 
=  (  T 
2.  / 
c  .1 
./  ( 
!>{K 
&(K 
-  X  X 
-XX 

1  )  = 
I)  = 

Alt 

Bl( 

Cl( 

X-5 

Y 

Z-S 

H*( 

H*( 

h*( 

.9* 

.9* 

.  9* 

100 

ISI 

ISI 

8  0 

C<  I 

0(  I 

=  (  [ 

2./ 

2./ 

PA* 

1*1 

If  1 

tit 

tit 

I)  = 

I)* 

I  =  S 

K  =  S 

K  =  S 


) 

»«  )  +  30( J  ) 

(YC(«)+ATAN2(2.*X*Y»T+X*X-Y+Y)/XK) 

+X*X-Y*Y)**2+(2.*X*Y)**2 

XK*X* (  T+  X*  X+ Y*Y  )/DEN 

XK*Y* ( T-X*  x-Y* Y ) /DEN 

XXTXX+YX+YX) 

)*COS( Yl  ) 

)  *  S I N  (  Y 1  ) 

*  X  /:  (  K, )  -  Y  X  *  Y  Z  (  H  ) 

*YZK  )+YX*XZ  (K) 

Aid)*  (jC(  I*1»K  )-S0(  I-l.K  )  ) 

C1(K)*(S0(I» K+1)-S0( I,K-1) ) 

I )*<G< I  +  1>J»K )-G( 1-1 tJ,K)  ) 

J)*(S(I»J-1»K)-G(I,J+1,K)) 

K) *{ G( I, J, K+l )-G( I, J»K-1) ) 

X (  I) *GY 


Z  (  1  )  *GY 
XX*U-YX 
(  P*YX+0 
(  P*AZ-C 
U1+.1*U 
V  1  ♦  .  1  *  V 
*  1  +  .  1**> 

K  =  K  <  ,  M 
GN (1»2* 

&N(1,2* 
1=2»NX 
) 

»K)+30( 
+X*X-Y* 
XK*X*  (  T 
XK*Y*( T 
( YC(K)+ 
)=G1(I  . 
)=G(  1,1 
K)=G(I  » 
K)=ZSYM 
Al( I) *< 
C1(K)*( 
0  (  I  + 1  #  K 
0( I»K+I 
0<  I+UK 


*  V) 

*AZ)*U+(P*XX+2*BZ)*V)+G*w 

*YX)*U+(P*SZ-0*XX)*V)+P*« 


(MZ-K)-l ) 
(  K-KK  )-l ) 


2  ) 

Y  )**? 

+  X*X  + 

-x*x- 

AT  AN2 
1) 
»  K) 
3>K  )- 
*G  (  I, 
S0(  1  + 
S0(  I. 
)-S0( 
SJ-SO 
+  IS)- 


+  (  2.*X*Y  >**2 
Y*Y  )/0EN 
Y*Y  )/DEM 
<2.*X*Y»T+X*X-Y*Y) /X<  ) 


(YX*Ul-XX*(Vl*C0S(Yl)-Wl*SIN(Yl))/Cij(<))/61(2) 

1>K  ) 

1,K)-S0(  1-1, K)  ) 

K+IS)-SO (I»K-IT) ) 

I*K)-30( I*K)+S0(I-1»K)+A3{ I )*OSI 

(I»K)-SO(I,K)+SO( It K-IT)+C3(K)*DSK 

S0(I-1  »<+!$)- SO  (H-1,K-[T)  +  S0(I-1,«-IT) 


SXX(I)=A<:(l)*liSlI 
SZ2(I)=C2(«)*DSKK 
S*Z{I)=T1*A1(1J*C1(K)*0SIK 

eo  cnMiNiit 

DQ     <*C      J  =  <:',<<Y 

X  =  A  U  (  2  ) 

Y  =  :•  0  (  2  ,  K  )  +  -i  0  (  J  ) 

Yl  =  ?A'MYC('<)  +  ATAN2(2.*x*Y,T  +  X*X-Y*Y)/Xi<) 

0tN=(T-K*y-Y'lcY)**?+(^.*X*T')**2 

>  <  =  t./>,Kn*(  f+x*x+y*y)/den 

YX=2./XX*Y*(T-X*X-Y*Y)/DEN 

P=C5(K)*CaS(Yl) 

C  =  C  I  (  K  )  *  M  N  (  Y  1  ) 

GY=Bl(J)*(G(2tJ-i»K)-G(2*J+i*K)) 

G?(1,J)=61(1»J) 

G 1 (1»J)-G(1»J»K) 

G(1,J»K)=GO»J>«)-<XX*01+YX/C4(K)*<?*V1-Q*W1)+SX(2)*GY)/AI{2) 

G(l,JtK)=Z5YM*b(l,J,K) 

NL=h  x-i 

X  *  A  0  (  N  X  ) 

Y=SO(NX,K)+B0(J) 

fl=PA*( YC( K) +U AN2(2.*X*Y,T+X*X-Y*Y )  /  X  K) 

CEN=(T+X*X-Y*Y)**2+(2.*X*Y)**2 

XXsg. /XK*X* (T+X*X+Y*Y ) /OtN 

YX=2./XK*Y*(1-X*X-Y*Y)/0EN 

P*C5(K)*C0S(Yi) 

Q  =  C5(K)*SIiN(Yl) 

GY*6HJ)*(G(NX,J-1,K)-G(NX,J  +  1,K)) 

G2(MX» J  )  =G1 ( MX, J  ) 

G1(MX,J) =G(MX, J,K) 

&(MX,  J,K)=G(Nl>  J>K)+(XX*U1+YX/C4(K)  *(  P*Vl-Q*  t*l )  + SX  (  NX  )  *t>  Y  ) /A  1  (  N  X  ) 

G( MX, J.<)  =  IS YM*G (MX, J ,K  ) 
9  0    CONTINUE 

G(1,1,<)=G(1,2,K)+G(2,1>K)-G(2,2,K)*ZSYM 

G(MX,i,K)=G(MX,2,K  )+G (NX , 1,< )-G ( NX , 2. K) *Z S YM 

E=G(1TE2(K),KY,K)-G(ITE1(K),KY,<) 

GU,KY  +  l,K)*G(MX,KY-l,K)-E*ZSYtf 

G(MX,KY+1,K)=G(1»KY-1,K)+E*ZSYM 

G(MX,KY»K)=G(1»KY»K)+E*ZSYM 

IF     ( iS*Ll .GT.O)     CALL     YSwEEP     (K) 
100    CTNTINLt 

IF     (KSYM.NE.O)     GG     TG     110 

JT  — 1 

G3    TG     20 
110    CONTINUE 

PR=1. 2*FK/ AAQ 

RM=1.2*rtM/NM/AA0 

GM=GM/NM 

RETURN 

END 
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2G 


SUb 
THE 
AND 
»HT 
COM 

1,  AG 

2,  CI 
32(1 
<tHA, 

COM 
1KG 

COM 
19), 
2,01 
L=K 
JL= 
IF 
C  (  1 
D(  I 
DO 
kj( 
Rl( 
Gl( 
Gl( 
CGN 
J  =  J 
13  = 
CON 
BC« 
DO 
AB  = 
AC  = 
X  =  A 
Y  =  S 
DEN 
XX  = 
YX  = 
FH  = 
H  =  R 
AZ  = 
8Z  = 
E=H 
E  =  H 
XXX 
YYY 
AA  = 
Yl  = 


F  0  U  I  1  N  t  Y  S  W 
E30ATIJNS 
SUPEKSJNIC 
ATED  DIFFE K 
MC'N  G(129,2 
(129)>A1< 12 
(17)  ,C2< 1 7) 
7  )  ,  Y  Z  Z  <  1 7  )  , 
CA,SA,FMACH 
MDN  /CtL/  P 

MJN  /bW?/  G 
S  Z  Z  ( 1  2  9  )  ,  R  0 

,  Q  2 ,  n  ^ 


(FWAC 
1-1)  = 
1-1)  = 
10  1  = 
I  )=1. 
I  )=1. 
I.  1)  = 
I ,  J  I  - 
T  INUt 
1 

12 

TINUfc" 
-Tl*rt 
fcG  1  = 
-Tl  +  A 
T1*A1 
G<  1) 
0(  I.K 
=  (T  +  X 
2./XK 
2./XK 
1  .-kg 
0(1)/ 
-  X  X  *  X 
-XX*  Y 
*AZ 
*BZ 
=  2.  /X 
=  <•.  /X 
P0(  I) 
YC(K) 


h.GE.l 

c. 

0. 

Il»  12 


G( 1,1, 
1)=G(I 


EEP  (K) 

FOP  G  AkE  SGLVEU  HbPE  »FOR  MIXED  SUBSONIC 

FLOW  ,EY  »0*  PELAXATIJN  ,AND  BY  USING  A 
ENCE  SCHctfF 

fc,17),SC(129,i7),E0(17),IV(129,17),ITEi(17),ITE2(17) 
cy)#A2(12^),A3(129),BC(26),^l(^t),82(2  6),6  3(2fc)»-'(17) 
,C3(17),C4(17),C5(17),XC(17),XZ(17),XZZ(17),YC(17)»Y 
KSYM,NX,NlY,NZ>KTEl,KrE2*ISYM,SCAL,SCALZ,XK,0r«EGA,ALP 
,T,PA, ABLAOE 
l»P2.P3,BtTA»FK»Fvl»DG,GM.[N,S,Ul,Vl»Ul.lK,jK,KR»io,JG, 

K129,26)»G?(12  9,26),SX(129),SZ(129),SXX(129),SXZ(l2 

(129),Rl(12^),C(129),D(129),Il,I?.L»»vlx,KY,MY,Tl,AAG 


. )  Jl=3 


L) 

f  Jl-1 


L  ) 


1( J  )*C1(K  ) 

11,13 

kd*b: 
(I  )*c 


k  j) 

(K  ) 


)  +  B0<  J  ) 

*X-Y*Y  )**2  +  ( 2.*X*Y )**2 

*X*( T+X*X  +  Y*Y  )/DEM 

*Y*( T-X*X-Y*Y)/DEN 

(I  )  +  XX*XX+ YX*YX 

FH 

Z(K)-YX*YZ(K) 

Z(< )+YX*XZ(K) 


3.*X*X+Y*Y)*DtN-4.*X*X*(T+X*X+Y*Y)**2)/DEN**2 
(DEN+2.*(T+X*X+Y*Y)*(T-X*X-Y*Y))/DEN**2 


K*{ (T+ 

K*X*Y* 
*Cbi K) 
+ATAN2(2.*X*Y,T+X*X-Y*Y)/XK 


-t  b- 


c 


?-X 

X«  Y 
) 

) 


l"       Tf-c      f-.'LL.    hJf' 
«  ?  I>  1  '      A  ■<  t     i)Eh 
RESPECT     1  [j    THE 

A  =  H*  X x 

f  =  h  *  y  a 

A<=h**.*((Y»v* 

A  Y  =  H  *  *  >  *  (  -  •?  .  *  >. 

P=AA*Cdi  I  Y  i  *  P  A 

C=AA*jIMYl*PA 

PX=-Q*YX 

PY=->J*X\ 

P7=-P*PA-w*F 

QX=P*PA*YX 

QY=P*?A*XX 

GZ=-0+P*FA*F 

E  X=- (  AX* ^Z (k  )  + 

bY=-(AY*X/(K)- 

EZ=-( A*x^Z ( K )+ 

FZ=-A*YZZ(<)+2 

YXJ=P*>3  +  Q*E 

Y  YJ  =  I-  n+^F 
7  x  j  =  P  *  E  -  0  *  <3 
ZYJ=P*F-J*A 
YXU  =  PX*!i  +  GX*E  + 
YXV=PY+ri+QY*E- 
YXW*PZ*d+QZ*E+ 

Y  ru=  h<  * A  +  ux  *  E  + 

YYv  =  PY*A  +  ;:r*F  + 
YYW*PZ*A+QZ*F+ 

Z  Xv  =  PY*t-G Y*^  + 
/X'«  =  PZ*E-QZ*P  + 
ZYU=  PX*F-3X*A- 
Z YV  =  P  Y*F- J Y* A  + 
ZYw=PZ*F-QZ* A+ 
DGI=G( 1+1, J, L) 
DGJ»G(I#J+1#L) 
DGK«G(I*J#L+1) 
DGII=G(I+l>JfL 
DGJJ=G(I*J+1»L 
UGKK=G(IfJ»L+l 
DG1 J  =  G( I+lf  J  +  l 
DGIK=G(I+1»J»L 
DGJK  =  G(  1  » J  +  1 »  L 
bX=  /  1  (  I  )  *UG  1 
GY=-31 (J  )*u^J 
GZ  =  C1U)*DGK 
UP  =  GX-SX;(I)*GY 
V«  =  GY 

WB=GZ-SZ(I )*GY 
U=UB+C A+xX+SA*YX* 
V=VB+S A*XX*P/C4 ( K 


LINES    iHf    COFFFiCItNTb    Of 
r:     >  A  S     a  E  L  L     AS     F  H  f  I  S     D  c  K  I  V  A  I 
*    CjC^OINATES 


x  *  *  <?  )  *  X  X  X  +  c  ,  * )  <  *  y  *  *  Y  Y  Y  ) 
X*XX>  +  (YX**2-XX**Z  )  *  Y  Y  Y ) 


THt    J 

IV  Pi 


A  C  "]  fi 

win 


A\ 


AY* 

AX* 
t*Y 

*xz 


P*A 
P*  A 
Q*fc 
F  *A 
P*  A 
0*F 
P*E 
P*E 
P*b 
P*b 
P*E 
P*F 
-G( 
-Gl 
-Gl 
)-G 
)-G 
)-G 
#  L  ) 
+  1  ) 
+  1  ) 


YZ(K) ) 

mo ) 

l  Z(K) ) 

Z  (K  ) 


+  0  *  E  X 

+  Q*EY 

-G*E  f 
+  Q  *  £  X 

-o*Ar 

+  Q*  A  X 

-0*Ax 
-  Q  *  A  Y 

-1 . J» L  ) 
ItJ-1) 
I  ,J  ) 

I»JfL)-G(IfJfL)+G(I-l»JfL) 
If Jf L)-G(1»J»L)+G(lf J-l.L) 
I, J  »L  )-G( It  J, I )+G( I. J,  L-l ) 
G(  1-1»  J  +  lt L )-G( 1+ 1, J-li L)  + 
G(I+1»J» L-l )-G( I-l# J. L+i )+ 
G(  I  .  J- 1 . L  +  1 )-&( It J  +  l » L-l  )♦ 


+  A2(  I 
-i3(  J 
+  C3(K 
G  (  1  -  l 
G  (  I—  1 
G  (  I  .  J 


)*0G1 

)  *L)GJ 
)  *0GK 
» J-i  t L ) 
»  J • L-  l  ) 
-1  »  L-l) 


P/C4U) 

)-CA*YX 


*  =  w 
UC  = 
vO  = 
bG  = 
QQ  = 

vo= 

1*0  = 
AA  = 

cc  = 

F  XX 
F  YY 
FZZ 
FXY 
FXZ 
FY2 
BU  = 
3V  = 
BW  = 
0U  = 
0V  = 
0'*'  = 

cu- 

CV  = 
1YW* 

cw= 

FXX 

FY* 

F  Y  f 

II)* 

FXY 
FYZ 
AV  = 
UU  = 
I  V  = 
UW  = 

vv  = 

v  W  = 
W'W  = 
M  X 
AZZ 
AXZ 
R  =  - 
1CV) 
AXT 
AY  I 

A/r 

A=* 
AXT 
AYT 
AZT 


b  +  C  A  *  X  I 
A*U-3*  V 
Y  X  J  *  u  +  Y 
ZX J* J+7 
u  C  *  *  <i + V 
VO-OMEG 
WG  +  'JffcG 
D  I  M  (  A  A  0 
U  0  *  U  0  +  V 
=  AUf  YX 
=  b  *  6  +  Y  Y 
=  C  *  w  +  F  * 
=  -&*  t$  +  Y 
=  0  *  Y  X  J  + 
=  (-*Y  Y  J  + 
A*UJ+>  X 
-B*Uo+ Y 
Q*VC+P* 
A*AX-b* 
-  A  *  A  Y  -  fa 

rx  j*wX  +■ 

UO* (AX* 

w) 

oO* ( -AY 
Fl>  ) 

VC  * (  j>* 
=  F  XX  +  k  1 

=  f y r  +  1 . 
=  SX(  H* 
SZ(1)*F 

=FXY-bZ 
=F YZ-bX 
EV-SX (I 

iiu*riu 
BU*AV 
bU*4w 

AV*AV 
Av*  ti* 
i?  W  *  b  W 
=F XX*t A 
=  F/LZ*AA 
=2. *(FX 
( AX A*bX 
*V3+( AA 
=  A  5  S  (  6  U 
=A3S( AV 
=  A  3  S  (  B  *l 
0(  1  )*bE 
=  A*A  XT 
=  A*A  Y  f 
=  A* AZT 


(K)  +  SA*(G  +  P*YZ(<))/CMK) 

r/J*\/  +  C*w 

YJ*V  +  P**i 

)**<>  +  wC**? 

A  *  P  /  C  4  (  K  ) 

A*  -/C^t  (K  ) 

>.2:*Ow-2.*Cj',1EGA*(3*U  +  A*V))) 

0 *vo+wo**c 

J*>xj+Zxj*ZXJ 

j  *  y  r  j  +  ?  y  vi  *  /  y  j 

p 

xj*yyj+zxj*2yj 

P*Z  Y  J 

J* VO+ZXJ*W0 

YJ*VO+ZYJ**0 

40 

AY+YXJ*YXU+YYJ*YX\/+3*YXV1+ZXJ*ZXU+7YJ*ZXV+F*2Xw 

*AX+YXj*YYG+YYJ*YYV+y*YYw+ZXj*ZYU+ZYJ*ZYV+P*ZY* 

fYJ*t.Y  +  C*GZ  +  ZXJ*PX+ZYJ*PY  +  P*PZ 

3U+AY*&V)+VO*(YXU*BU+YXV*3V+YXW*tf*)+WQ*(ZXU*3U+ZXV*3V+ZX 

*ts  L:  +  AXtbV)+Va*(Yfu*BU+YYV*3y/  +  YY"w*oW)+wC*(ZYL*DL+ZYV*bV  +  Z 


3  U  +  Q  Y*  6  V  +  3  Z  *  a  W  )  +  WG*  (  ?  X*  6U  +  P  Y*  t)  V+  P  Z  *  B  V ) 
(  I  )  -  k  0  (  I  ) 
-»0(  1  J 

X?*rXX+FYY+SZ(i)**2*FZZ-2.*iX(I)*FXY-2.*<-Z(I)*FY 

X.Z 

(  I  )  *  FXZ-SX (  1  )  *FXX 

(  1 ) *FX7-SZ (  1  )*FZZ 

)*BU-S2 ( 1 )*8w 


7  +  /•'.*  S  X  ( 


-bU 

-«  * 

Z* Al-UW  ) 

<(I)+AZZ*S72(l)+AXZ*3XZ(I))*GY+TlM(AA*Ji-CU*ct) 

*DW-Cw)*WB) 

*  A  1  (  I  )  ) 
*»  L (J  >  ) 

*  C  1  (  K  )  ) 
rA*AA/AMAXl(A<T,AYT»AZT*l.-SO(I)) 


♦  (  a  A  *  D  V  - 


3  0 


40 


DC 


SUPERS 

IF  (0J 

a  y  *  =  i  k 

«f  Y  =  (  F 
A  Z  Z  *  A  Z 

AX Y=  (  f 
AX  Z  =A  K 
A  Y  7  =  (  F 
t  P  =  A  X  X. 
B  vi  =  A  X  X 
B*- A  XX 
R  =  A  X  X  * 
GG  TO 
CONT  IN 
N  S  =  N;>  + 

s  i  =  s  i  g 

IM=  l-i 

IMM»IM 

AX  X  =  UU 
AYY  =  M\J 
A  Z  Z  =  '*  w 
AX  f  =  b. 
AXZ  =  8. 
AYZ  =8. 
dx a  =  (F 
BYY=  (  F 
bZZ=(F 
BXY=(F 
bXZ=(F 
3YZ=(F 
AQ= AA/ 
DEL  TAG 
DGI  I  =  G 
DGJ J=G 
DGKK=G 
D  G I  J  =  G 
DGIK=G 
DGJK=G 
G  S  S  «  A  X 
B».5*< 
PP=  ACi* 
BM= A3* 
b=-AQ* 
1  ) 
R= ( AQ- 
CONTIN 
IF     (A3 

IP-  I 

J9  =  J 
K«!  =  K 
CGNTIN 


I-nJC     *3 1  NTS     AMD    jU^IMC    POINTS     A- 

. G  t . A  A )     GO     Tu    31 

x*A3  (  1  ) 

YY vH a-v V  )  *B2  (J  ) 

7*CZ ( K  ) 

XY*AA-U\/  )*(  Ae  +  AF) 

i  *  AC 

YZ*AA-VW)*(8C+3C) 


StPA^.ATFL 


-AXX 

D&II 

40 

UE 

1 

Ml. 

Flx( 

-IFI 

*A^( 
*B2< 

*C2  ( 
*SI* 

ax  *u 
T  Y  *U 
ZZ*Q 
XY*J 
XZ*G 
YZ*0 

au 

=  3Xx 

(  i>j 
( I  ,j 
(  i,  j 
(i#  j 
(  i»j 
(It  j 

x  *  b  G 
Ag-1 
BAX- 
3  XX- 
(  3XX 

1.)* 

Lit 
S(R) 


ut 


-Q1*(4YY+AZZ) 

+AYY*0GJ  J  +  AZZ*OG><N+AXY*QGI  j  +  A  YZ  *C&  JK  + A  X  7*  OG  I<+ R 


.1)  ) 
i  1  ) 
X(  SI  ) 

1 ) 

J) 

K) 

JV*  A  E 

J.v*AC 

3C 

J-UU  )  *A2 ( I ) 

3-VV  )*S2( J) 

J-nW  )*C2  (K  ) 

J-UV  )* ( AB+ AF ) 

Q-0*  )  *( AC  +  AC  ) 

i-M'n  )*  (BC  +  bC  ) 

*DGII+BYY*DGJJ+8ZZ*DGKK+BXY*DGIJ+8YZ*DGJK+BXZ*DGIK 

»L  )-G(XM*J*l)-G(  IM> J»L)+G(IMM*J«L)  +  A3( I)*DGI 

,L)-G(I,J-l»L)-G(i,J-l,L)+Gl(l,J-2)-P3(J)*l&J 

»U-G(  I»  JjL-D-GCI*  J»L-1)+G2{  I»  J)+C3(K  )*DGK 

»L)-G(IM»J»L)-G(IfJ-l,L)  +  G(n»J-l,L) 

iL)-G(If JjL-l)-G(IMf J,L)+G(IM,J,L-1) 

.  L)-G(l»J,L-l)-G(l,J-I,U+G(I,J-i,L-l) 

II+AYY*0GJJ+AZZ*0GKK+AXY*0GIJ+AYZ*DGJK+AXZ*D3IK 

.  ) *(  AXX  +  AXX  +  AXY  +  AXZ  ) 

(l.-SI  )  *3 

(l.+SI)*B 

+BXX+Q2* (BYY+3ZZ ) )+( AQ-1. )*( 2.*( AXX+AYY+AZZJ+AXY+AYZ+AXZ 

GSS+AQ*DELTAG+ft 

.LE.  £BS<  FP  )  )  GO  TO  50 


-8  6- 


6C 


7C 


to 


QC 


R  =  R-AYT*(GK  It  J-1)-G(I»  J-lf  L)  )-AZT*(Gl<  if  J  )-G ( I f J » L-l ) ) 

Rtf  =  R!*  +  A8$(R) 

N  M  =  N  i"l  + 1 

b=B  -AAT-AfT-AZT 

BM»BM+A<T 

B  =  l . / < tf-6*!*C ( 1-1)  ) 

C(  1  )=d*3P 

D(  I  )  =  B*(K-BI"i*0  (  1-1  )  ) 

CONTINUE 

CG  =  G. 

1  =  13 

DO  80  M  =  1 1 # I  3 

CG«D(1 )-C(I)*CG 

GM=Gtf+Ai3S (CG  ) 

IF  (A3S(CG).Lt.A6S(DG))  GO  TO  7  0 

DG=CG 

IG=I 

JG=  J 

KG  =  K 

CONTINUE 

G2(I»J)=G1( I* J) 

Gl(  If  J  )  =  G(  I, J»L) 

G( I» J»  L)=G( If  J  >L  )-CG 

1=1-1 

CONTINUE 

J  =  J  +  1 

IF  (J-KY)  20»9CfllO 

con  riNut 

IF  (I2.GT 


I3=ITE2(K) 
13  =  LX 


ITE2  (K)  ) 

IF  (  ITF2(K)  .EQ.HX) 

CO  100  I  =  I  i  »  13 

LV=  IA6b ( 1-IA3S (  I  V  (  If K  )  )  ) 

kO  (  I )  =  A  M  I  NO  (  L  V  »  I  A  B  S  (  I  V  (  I »  K  )  )  ) 

Pl(  I)=LV 
IOC  CONTINUE 

GO  TO  20 
110  CONTINUE 

I=LX+i 

IO=NX+2-I3 

IF  (K.GT.KTE2)  GC  TO  130 

DO  120  I =1  Of  I  3 

X=AO(  I) 

Y=SO(IfK) 

DEN=(T+X*X-Y*Y)**2+(2.*X*Y)**2 

XX*2./XK*X*(T+X*X+Y*Y)/DEN 

YX=2./XK*Y*(T-X*X-Y*Y)/DEN 

A=l .-R0(  I)  +XX*XX+YX*YX 

H»RO(I )/A 

A7=-XX*X2(<)-YX*YZ(K) 

BZ=-XX*YZ( K) +YX+XZ (K) 

A=RO( I )*Cb(K) 

Yl*YC<K)+ATAN2(2.*X*YfT+X*X-Y*Y)/XK 


-8<i- 


120 


130 


14C 


p  =  i 
C=A 
FZZ 

r  y  r 

hX  < 

f  x  y 

FXZ 

fyz 

AX  = 

AY  = 

AZ  = 

BY  = 

DG1 

DGK 

V  =  S 

U  =  A 

*  =  C 

6(1 

CGN 

1  =  1 

IF 

E  =G 

£C( 

CON 

IF 

1  =  1 

t  =  0 

IF 

b  =  E 

CON 

M*N 

6(1 

G  (  M 

G2( 

Gl( 

G(M 

GO 

END 


+ c  a  s  ( 

*  s  I  N  ( 

=  A*  A 
=  H*r^* 
=  H  ^H* 
=  H*H* 
=  H*FZ 
=  H  *  F  Z 
Fx  f-i 
F  Y  Y-j 
F  f  Z  -  .s 
AY-iX 
=  G  (  I  + 
=  G(  I  . 
A*X  A* 
1(1)* 
1(K)* 

>ky  +  ; 

t  i  n  u  e 
c 

(  1 0  .  N 

(13>< 

L)  =  eo 

I  INbE 
(I. It 
-1 


^*  Yi  ) 

P4*y 1) 

( Y  <* Y A 
I  X  X  *  A  X 
(-XX* Y 
/  *  A  7 
I*  it 
x (  I  )*F 
x( 1 )*F 

A ( I ) *F 

(  I  )  v-  A  X 
i»<Y,L 
KY» L  +1 
P/C^  (K 
Cbi+C A 
DGK+CA 
.  L  )  =  G  ( 


+FZZ*(8Z*BZ+XX*XX))+1.-K0(I) 
+FZZ*lAZ*AZ+YX*YX))+i.-fiG(I) 

X  +  FZZ*(AZ=MZ+AX*YX)  ) 


X  X  -  >  Z  (  I  )  *  F  X  1 

XY-S7(  I  )*FY7 

XZ-SZ  (  I  )*i-ZZ 

-  S  Z  (  I  )  *  A I 

)-C(  I-1#KY»L  ) 

)-G2  (  I, k Y) 

)-CA*YX-JMEGA*XX/C^(K  ) 

*XX  +  SA*YX*P/C<»  (K  )-Q>lEGA*YX/C4(K) 

*xZ(K)  +  $A*(Q  +  P*YZ(K))/CMK)-DPEGA*YZU)/CMK) 

I»KY-ltL)+(AX*U+AY*V+AZ*W)/(PY*Hl(KY)) 


E.ITE1 (K)  )     GO    TO 

Y>  L )-G (  IO»KY, L  ) 
(L )+P3*( t-£0( L ) ) 

.11)     ^FTUKN 


1  3i 


(IVli*K).NE.l)    GO    TC    l<tO 

G(L  ) 

TINUt 

X+2-I 

,KY  +  i,L  ) =G(  l",KY-l»  L  )-E 

>  KY  +  W  L ) =G(  I ,K Y-lj  L )  +  b 

M>KY)*GL (M»KY) 

M,KY)=G(M>KY>1) 

>KY.L)=G(1*KY>L)+E 

TO    130 


SUBROUTINE     VELO     (K,L»SV,3*>CP>  X,Y) 

CALCULATES    SURFACE     VELOCITY    AND    PRESSURE    COEFFICIENT 

C  OIMON     G 1 129, 2  6, 17 ) j  S C ( 129, 1 7 ) >  EO ( 1 7 ) » I  V ( 1 2 9  » 1? ) > I T E  1  ( 1 7 ) j I T  F 2 ( 1 7 ) 

l>A0(12n,Al(i2^)»A2(12<5)»A3(12S<),BC(26)^i(2b),a2(2c).^3(2b)»Z(17) 
2.CK17),C2(i7),C3(17),Ct(17),C5(17),xC(17),xZ(17),xZZ(l7),rC(i7),Y 

3Z(17)»YZZ(17)»KSYM>NX»iNiY#NZ*KTEl»KTE2»ISYi,'»bCAL»SCAL?»XK»J^tGi.ALr 


-90- 


h  H  4  »  C  A  >  ;WI-  1aCh»T»PA>ABLA0E 
UIMENSIUN     S  v  (  1  )  >     S  M ( 1 ) »     CP(x)»     <  (  1  )  » 
110=  I  T  t  1  (  k  ) 
I  2  0  =  1 T 1 2  ( K  ) 
J  =  N  Y  +  1 

01=  .  ? *  Frl  A C  -1**2 
T 1  =  1 .  /  (  .  ?  +  f- fc.  i  C  H  *  *  2  ) 
B=C5 U  ) 

C  J     1  C     1  =  1 1  J >  I  £  0 
A=AG( I ) 
Q  =  S  G  (  I  J  K  ) 

DEN=(T+A**2-Q**2)**2+(2.*A*G)**2 
XV=?,/XK*A*( T+A**2+n**2)/DEN 
YX=2./XK*Q* ( T-A**2-0**2) /DEN 

CL.-V    v*V    *    -A.  V  V  £  V  \ 


Y(  1) 


YX=2./X>, 

Fh=  XX*X  UYU  YX 

H  =  0. 

IF     (IV(I*K).NE.O)     H  =  1 .  /  F  H 

AZ 

B 

D 

1>X=  Al(  i )  *0S  1 

SZ  =  C1  (<  )*DS«- 

Y1  =  YC(K)+AU,>I  2(2.*  A*  G»T  +  A*A-0*0)/XK 

P=B*CGi ( PA*Y1  ) 

C  =  3*SIN ( PA*Yi  ) 

Of 


Z  =  -x  K*x  ?(K)-Y  X*YZ (K ) 
iZ  =  -Xa*YZ(K)+YX*XZ(K) 
>SI*SO(  I  +  1»K  >- SO ( 1-1»K  ) 

SK=SO(l*K+i)-SC(I#K-l) 

k  s  a  i  r  n  *  n  ^  i 


A    V    L    /    -   A  \„    I   l\   #    '    n  L    _»  \J    %    J 

Y  (  I  )  =  Y  C  (  K  )  ♦  A  T  A  N  2  ( 
A2( I )=XZ (K )-YX*SZ 
A3(  I )=YZ  IK )  +  XX*SZ 


10  CONTINUE 
RETUPN 
ENO 
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1C 


20 


30 
40 


20 


sua 

PLG 

JIM 

PA  T 
1WR 
».PI 

DO 
LIN 

FDE 
A^A 
CMM 
CPM 
DO 
A  1A 
C  M  M 
COM 
CM  A 

DC 

XP<? 

Kl  = 

Kl  = 

LIN 

K?  = 

K2  = 

LIN 

K3  = 

K3  = 

IP 

K3  = 

GO 

LIN 

K4* 

LIN 

JJ  = 

W*I 

LIN 

L  IN 

LIN 

LIN 

CJN 

RET 


TS 
ENS 

A      K. 

1T  = 

Tr 
lj 

E(  I 
N  =  l 
X  =  ) 

Ax  = 

AK  = 

2D 
X  =  A 
AX  = 
a  K  = 
X  =  C 
50 
AC  = 
(  5  S» 
Mn 

E  (K 

41 

KIN 

E  (K 

(  3o 

MIN 

IK3 

1. 

TG 

E  (K 

1 

E  (< 

0 

Tfc 

E(K 

E  (K 

E  (K 

E  (* 

TIN 

URN 


I  I  n  t    C  P  L  0  T     (I1»12*X,},D*FMACH) 

CP    aT    EjUAL     ir:TE-VAL>     IN     ME     "IAPPE 

ION    <JDc(3)j     LINE ( 100) t     X(l),     3(1) 

iluc  /  1h     ,  3ri+,lHQ/ 

c 

( I  k  K  I T  »  b  0  ) 

1=1* 100 

) =KJOt  (  1  ) 

.  /  b  (  I  P  ) 


[=11.12 

«  A  X  I  (  A  ,"i  A  X  ,  A  8  S  (  X  (  1  )  )  ) 

A.m]NKC*.MAXjD(I)) 

AMAXKCPMAX.D(I)) 

PMAX-CMM  A  X 

1=11*12 

E  l)  E  N  *  r<  (   I    ) 

./AfUX)*A9S(X(I))-K»l. 

0  (  K  1  t  1  0  J  ) 

1  )«KQDE( 2) 

J(K2. 100  ) 

1  )=<JOc  (3  ) 

./CNUX)*(D(I)-CM--1AX)+l 

0(K3»40) 

.oE.l)     GO    TO    30 

3  )=KODt (2) 


4)=«JDE(  3  ) 

(  IwPIT»70) 
1)=K3DE(1) 
2)=KQDfc( 1) 

3)=K00E( 1 ) 
h  )  =KQDE ( I  ) 
UE 


E  L  A  N  F 
J(  J   ) 


XFk AC,  X (  I ), JJ, LINE 


fO 
70 


FOR 
FOR 

ENO 


MAT     (  lx»  It  3  X»  3HX/C  »4Xj  3-liviC  GMP  »  ^  K  .  5HM0SGN*  3H    QN) 
MAT     ( IX , Fb . 2»F9.  3, I3>2X,iCQAl ) 
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THE 

PCI 

COM 

1>A0 

2, CI 

3Z(  1 

<.HA» 

DM 

Cl  = 

A  =  C 

11  = 

12  = 
K  Y  = 
wRI 
wRI 
DG 
DSI 
DSK 
bX  = 

sz  = 

DO 
X  =  A 

Y  =  S 
DEN 
XX  = 
Yx  = 
FH  = 
H  =  0 
IF 
AZ  = 
BZ  = 
Yl  = 
P  =  A 
0=A 
DGI 
DGJ 
DGK 
GX  = 
GY  = 
U  =  G 

V  =  G 
W  =  C 
UO« 
VC  = 

wo= 

Q  =  U 


UT  I  ^t 
PttC 

,  TH 
N  G(l 
2s>)  ,A 
7)  ,C2 
,  YZ2  ( 
,  SA,F 
SIGN 
/f-MAC 
K) 


i  p  e  e  d 

AND  T  H 
EN  THE 
29, ^6» 
1( 12$) 
( 17) ,C 
17), KS 
>l  A  C  H  ♦  T 
IS  (24) 
H  *  *  2  ♦  0 


(K  ) 
t       MACH    rtUMBEP     ARE     CUMPUTFO     AT    EACH    GRID 

"ACH  MU^BERS  ARE  WRITTEN  ON  THF  01  T  F-UT 
17),  3  0(129,17)*  E0(17), IV(129,17)»ITEl(17),irE2(  17) 
»A2(129),A3(129),B0(26),B1(25),B2(26), 83(26), Z  ( 1  7  ) 
3(17),Ct(17)»C5(17),XC(17),XZ(17),XZZU7),YC(l7),Y 
YN.NX»NY.NZ,KTEl,KTE2,ISrM,SCAL,SCALZ»X<,G'/EC;^,ALHJ 
, PA, ARLADE 


.2 


U 
(6*30) 

(b,<tO)     Ff.ACH 

1=11,  12 
3(I+1»K)-S0(I-1,K) 
0(1, K  +  l  )-SC(  I,K-1) 
(  1)*D3I 
(K  )*OSK 

J  =  <-,<Y 
I) 

i,K)  +  ao( j  ) 

T+X*X-Y*Y )**?+( 2. *X*Y)**2 
/A**X*(T+>*X+Y*Y)/DEN 
/XK*Y* ( T-X*X-Y*Y  )/OEN 
**<^+Y<**2 


(FH.GT. .16-09)    H=1./FH 

-XX*XZ(K)-YX*YZ (K  ) 

-X  X*YZ( K)  +f  X*XZ ( K) 

YC(K)+ATAN2(2.*X*Y,T+X*X-Y*Y)/XK 

*C0S(PA*Y1 ) 

*SIN( FA*Y1) 

=  G(  1+1,J,K)-G(  1-1 »  J,<) 

=  G(  I, J+ 1,K )-G( I, J-l »K ) 

=  G(  I» J, K+l )-G ( I, J, K-l ) 

Al  (  I)*DGI 

-BK  J  )*0GJ 

X-SX*GY+CA*XX+SA*YX*P/C4(K) 

Y+SA*XX*P/C<t  (K  )-C A*YX 

1(K)*QG*-SZ*GY  +  CA*)<ZU>+:>A=MC+P*YZ(K))/C4('<.) 

H*( XX*U-YX*V  ) 

H*((P*YX+Q*AZ)*U+(R*XX+Q*BZ)*V)+0*W 

H*((P*AZ-0*YX)*U+(P*BZ-Q*XX)*V)+P*w 

0*UG  +  V0*\/0+W0*W0 


-93- 


IS(J)  =  S|jPT(U/(Q1-0.2*(Q-2.*0MEGA*HMYX*u+XX*V))))*]0CC. 
10  CONTINuc 

WRITE  (6*50)  (  IS ( J), J=2»KY) 
2  0  CONTINUE 

RETURN 


30  FORMAT  (1H1) 

40     FORMAT  (1SH     PRINTOUT     OF     SPEED, /,13H 

50    FORMAT  (2415  ) 
END 


FMACH 


4.  2,/) 


1C 


SUbR'JUTINE  FOPCF  (I1,I2,X,Y,CP,AL,  CHORD,  xr.,CL,CD,Ct'»Xt<,YK,Zl,PA) 

CALCULATES  SECTION  FORCE  COEFFICIENTS 

DIMENSION  X(l),  Y(l),  XK(1),  YK(1),  CP(1) 

RAD  =  57. 29578 

ALPHA= AL/R AD 

CL*0. 

CD  =  0. 

CM=0. 

N=I2-1 

DO  1C  1=11, N 

DX=(  X(  I  +  D-X  (I  )  WCHORD 

OY=(Y(  1  +  1 )- Y ( I  )  )  /CHORD/ZI 

XA=(.5*(X(I+1)+X(I))-XM)/CHQRD 

YA=.5*( Y(I+1 )+Y(I ) ) 

CPA«.5*(CP(I+1 )+CP ( I ) ) 

DX*DX*COS(PA*YA)+XK(I)*SIN(PA*YA)+DY 

DY«DY*<l.+YK(I)*SIN(PA*rA)**2/ZI) 

YA*((1.-PA)*YA+SIN(PA*YA)/ZI)/ CHORD 

DCL=-CPA*DX 

DCD=CPA*DY 

CL=CL+DCL 

CQ*CD+DCD 

C*1  =  CM  +  DCD*YA-DCL*XA 

DCL«CL*COS(ALPHA)-CD*SIN<ALPHA) 

CD=CL*SIN(ALPHA)+CD*COS(ALPHA) 

CL=DCL 

RETURN 

END 
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10 


sua 
ins, 

CAl 
DIM 
SPA 
CL  = 
CD  = 
CMP 
CM* 
CNY 
S  =  0 
N  =  K 
DO 
DZ  = 
AZ  = 
DZ  = 
AZ* 
CL  = 
CD  = 
CMP 
1D(K 
CMR 
CMY 
S=S 
A3L 
CL  = 
CD- 
CMP 
CMR 
CMY 
RET 
END 


iRJUFINE  TOT  FOR  (KTtl,KTE■i»CHLl*D^CL»i:0#iCh.Z>^C.C•J»lL»C•L,  >CMP,C 

CMY.»  PA>  ArtL  4Cf  ) 

CULATtS  TOTAL  FO^Ct  COEFFICIENTS 

ENSIGN  C  H  0  8 C ( 1  )  t     S  C  L (  1 ) »  S  C  0 ( 1 ) »     S  C 1 ( I  )  »  Z  ( 1  >  .  X  C ( i ) »  C b ( 1 ) 

N=Z(KTE?)-Z(KTE1J 

C. 

0. 

=  w  t 

=  0. 

=  0. 


TE2 
10 
.5* 
.5* 
DZ  + 
AZ  + 
CL  + 
CD  + 
=  CM 
)*( 
=  CM 

=  cm 

+  DZ 
ADF. 
CL/ 
CD/ 
=  CM 
=  (C 
»(C 
URN 


KTE1»N 

,/Cb(K+l)-l./C?(K) ) 

,/C5U  +  l)+l./C5(K) ) 

.-PA )* (Z(K+1)-Z( K) )/2. 

.-PA  )* (Z(K  +  1)+Z(K)  ) 12. 

*(SCL(K  +  1)*CH0kD(K+1 )  +  SCL(K)*CHORD{K)  ) 

*(SCD( K+l) *CHCRD(K  +  i)+SCO(K )*CHORD(K )  ) 

0Z*(C'H0RD(<  +  1)*(SCM(K  +  1)*CH0RD(K+1)-SCL(K+1)*XC(K*1))+CMi 

M(K)*CHORD(K)-SCL(K)*XC(K))) 

AZ*DZ*(SCL(K+1)*CHGkD(K+1)+SCL(K)*CHURD(K)) 

AZ*DZ*(SCD(<+l)*CHOkO(K+l)+SCO(K)*CHJHO(K)) 

CHJRD(K+1 )+CHORD(K ) ) 


SPAN/S**2 
+CMR) / (S*SPAN) 
+C^Y)/(S+SPAN) 


SU6R0UT 

HALVES 

COMMON 

1»A0(129 

2»C1(17) 

3Z(17), Y 

4HA,CA, S 

MX*NX+1 

KY=NY+1 

fY=NY+2 

MZ=NZ+3 

MXO=NX/ 


INE  PEF  IN 

MESH  SIZE 

&(129»26»17)*SC(129»l7)»EO(17)fIV(129,17)»ncl(17)*ITt2(17) 

),Al(129)>A2(129),A3(129),dQ(26>,BH26)»!'?(26)»63(26)»Z(l7) 

,C2(17)»C3(17)#C4{17)*C5(17)*XC(17)>XZ(17)jXZZ(17)jYC<17)#Y 

ZZ(17)»KSYM,NX,NY,NZ*KTEl*KTE2iISYM,SCAL»SCALZjXK,QMFc-.A,ALP 

A,FMACH,T,PA, AELADE 


2  +  1 


-9b- 


1C 


20 


30 


<i0 

50 

eo 


70 


80 


9C 


100 


l*.Z0 

KK  = 

DD 
J=N 
JJ  = 
I«M 
11  = 
G(  I 

I  =  1 

II  = 
IF 

J  =  J 
JJ  = 
If 

DO 

DO 

G<  I 

DO 

DC; 

G(I 

G(I 

CON 

MZM 

MZS 

COM 

DO 

DC 

G(  1 

IF 

MZS 

DO 

DO 

G(I 

NZM 

IF( 

MZS 

GO 

CON 

KK  = 

DO 

I  =  M 
IF 

II  = 
12- 
DO 
DSI 
DSK 
SX  = 
sz» 

DGI 
DGK 
R  =  A 


=  NZ  /2+3 
i 

6 

r 
K 

X 

I 


G     K  =  *  K  •  h  Z  j 
/2  +  1 

y 
o 

X 

»  J  J  »  K  )  =  G  (  1  »  J  t  K  ) 

1 

1-2 

l.GT  .0)     Gj    TL    20 

1 

J-2 

J.GT.G)     GO    TO    10 

C    J  =  1>KY>2 

0     I  =  ^ »  N  x  ,  2 

J.K)=.5*(G(I+l»J»K)+b(I-l»J.K)) 

J     I-ltMX 

0    J«2»NY»2 

J»K)=.5*(G(I»J+1»K)+G(I,J-1»K)) 

m  y  *  k  )  -  o . 

INUE 

MZO 

*MZ 

INUE 

o  j  =  i,n 

0    I=1#MX 

J,MZST )=G( I»J#MZK) 

MZST.tO.l)     GG    TO    100 

=  M  Z  S 1  - 1 

0    J  =  ]  *  M  Y 

0    1  =  1*  MX 

J»MZST)*0.5*(G(I#  Jt  MZ«)+G(  1»J>MZM-D) 

MZM-1 

ZST.fcQ.l )     GO    TO    iCO 

=MZST-1 

0    70 

INUE 


5  0    K=KK,MZ 

0+1 

K.LT.KTE1.0R.K.GT.KTE2) 

TEl(K) 

TE2(K) 

10    1  =  11*  I? 

S0( I+1»K)-S0( I-1#K ) 

S0( I><+1 )-S0( I»K-1) 

1(1) *OSI 

1(K)*DSK 

G(  I+1*KY»K )-G(  I-1»KY,K  ) 

G(  If KY»K  +  1 )-G(  I*KY#K-1) 

INOUjIVUjK)) 


,3    TO    12G 


-9o- 


110 


120 


13C 


140 

150 


AC  (  I 
SCC  I 
N=(T 
=  1.  I 
-2.1 
1.0- 
S/A 

=  -xx 
=  -xx 

P*Ct> 
=  YC( 

A*Cl 
A*SI 
Z  =  A* 
Y  =  H* 
X=H* 
Y=H* 
Z=H* 

=  FXY 

=  FYY 
=  FY? 
=  AY- 
Al(  I 
CKK 
S  A*X 
I»  <Y 
(  K)  = 
11 
1-1 
0. 

(IV 
EG(K 
NX  +  ? 
I»KY 
M,KY 

(  IV 
I»KY 

(  IV 
,K-1 
M»KY 
I»KY 
M»KY 

(  I. 
NT1N 
TURM 
D 


) 

fx.) 

+  X  *  X 

x**X 

ft  +  a  x 
*</( 

*YZ  ( 
(K  ) 
*  )  +  A 
l(  PA 
N(  PA 
A 

h*(  r 

H*  (X 

H*  (- 
FZZ* 
F  ZZ* 
-SX* 

-sx* 

-sx* 

jX*  A 
)*0G 
)*u(. 

A*  P/ 

f  liK 
&{  12 


(I»K 
) 

-I 

+  1*K 
+  1»K 

( i»K 

,K)  = 
(  I,K 
)  ) 

,K)  = 
-1*K 
-1»K 
GT.2 
UE 


-Y*Y)**2+(2.*X*Y)**2 

MT>X*X+Y*Y)/DEN 

*(T-X*X-Y*Y)/DtN 

*  x  x+  Y  x  * Y  x 

<)-YX*YZ(K) 
K)+/X*XZ(K) 


r  an?  ( 

*Y1) 

*Yi  ) 

X  *Y*  +P 
A*XX+F 

xx*  y  m 

AZ 
3Z 

FXX-SZ 
FX Y-SZ 
FX7-SZ 
X-SZ*A 
I+C A*X 
<+C A*X 
C4U)- 

>=G(    I. 

»K  Y»K  ) 


.*X*f, T+X*  X-Y* Y  )/ X' 


2  Z  *(  bZ**Z+XA*XX  )  )  +  1  .-K 
Z  Z*  (  AZ*AZ+Ya*  YX  >  )  +1  .-(5 
F  Z  Z*  (  AZ*BZ  +  AX*YX  )  ) 


*FYZ 

*  F  ZZ 

Z 

X+SA*YX*P/C4(K)-QMFGA*YX/C4(K) 

Z(K)  +  SA*(Q  +  P*YZ(K))/CMK)-0MEGA*YZ(K)/C4{'<) 

CA*YX-QMEGA*XX/C4 (  * ) 

KY-1>K)+(AX*U+AY*V+AZ*W)/(BY*61(KY)) 

-G  (  II  .*Y,K ) 


I.Nc.l)    GG    TJ    130 


) =  &( M.KY-1 . K )-£ 

)=G( 1 »KY-1 »K ) +t 

I.NE.-l)    GO    FO    14o 

.•>*G(I,KY,K-1)  +  .2:>*(&(I>KY.K+1)+G(M,KY,K  +  1>) 

+  D.LT.1)     G(I»KY»K)».5*G(I*KY*<  +  l)  +  .25*(G(I*KY*K-l')  +  G(Mf 

G(  I>  KY,K) 

)=.;*(G(I»KY,K,)+G(  I.KY-2,  K)  ) 

)«.5* (G(M»KYtK)+G( M»KY-?»K) ) 

)     GJ     Tu    120 
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Sl 

LV 

CI 

If 

IN 
K  = 
IF 

1C  K  = 
1  = 
J  = 
DS 
D  = 
IF 

2C  UF 
IF 

3C  U  = 
V  = 
GG 

40  U= 
v  =  ' 
GG 

bC  U  =  - 
V- 
GO 

60  1*. 
J*, 
DS 
IF 

70  DF: 

3 


R  C  L)  T  I N  c  5  P  L  I  F  (M»N,3»F*F?*FPP*FP'»P*i<M*\/M»KN*\/N»lw 
IC  SPLINE  FIT  wITh  ^ESCRIBED  END  CONDITIONS 
ENSIJN  id),  F(l)»  FP(1),  FP»(1),  Ff>PP(i) 

(  1 ' J  0  .  t  U  .  0  )  G  0  T  G  1  e  0 

MS  IN-M  ) 

(K-l)  13G,ie0*10 

M-M)  /  r< 

+  K 

S(J )-S(I ) 

IS 

(DS)  20*180*20 

(F(J )-F(I ) ) /DS 

(  K .  M  -  2  )  30*40,50 

5 
i.*{UF-V|Vl)/DS 

TO  dO 


I  N  0  ) 


OJ 


M 
I! 
1. 

DS*V.1 
TO  6  0 


U" 

v  = 

80  FP 
FP 

0  = 

v  = 

IF 

90  IF 

ICO  v  = 

GO 
110  v  = 

GG 
120  V» 
130  B- 

D  = 
140  DS 

U  = 


J  )- 

*bi 
(  J) 
(DS 

S 

b.* 
=  U 
)  =  V 
-0) 

DF  + 
-is) 
N-2 
*  VN 
U 


S(  I  ) 

)  1^0,190,70 

-F (  I  )  )  /  CS 

+DS+U) 

DF-V) 


*D5 
.*DF+DS*V 
(J-N)  63*90*60 
(KN-2)  100*110*120 
6.*VN-V) /U 
TO  13C 
N 

TG  13c 
DS*VN+F.JP(  I  )  )/  (l.+FP(  I)  ) 


S( J  )-S(  I  ) 
PP(I )-FP( I ) *v 
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FPPP( I )= (V-U)/ DS 

F  P  P  (  I )  ■  U 

FP{  I)>(F(J)-F(  I  )  )/DS-D5*(V+U+U)/6. 

v  =  u 
J-I 

I-I-K 

if    ( j-n >    1^0.1  50>i<»o 

150    I  =  N-K 

FPPP(N)*FPPP(I) 

FPP  (N  )  =  B 

FP(N)=Di-+u*(FPP(I)  +  3+B)/b. 

IND-1 

IF     (MODE)     130>180>160 
160    FPPP(J)  =  rQ:1 

V  =  FPP( J  ) 
170    I=J 

J  =J  +K 

DS=S( J)-S( I ) 

U=FPP( J  ) 

FPPP(J)=FPPP(I)  +  .t>*DS*(F(I)  +  F(J)-DS*DS*(J  +  V>/12.) 

V  =  U 

IF  (J-N)  17J»180»170 
180  CONTINUE 
RETURN 
END 


10 


dO 


30 
<tO 


50 


SUBROUTINE  INTPL  (MI »N I#SI t F I #M*N*S#F* FP# FPP# FPPP* MODE  1 

INTEkPClATIJN  OF  CUBIC  SPLINE  BY  TAYLCR  SEPIES 

DIMENSION  SHI),  FI(l),  S(l>,  F(l),  FP(1).  FPP(l),  i-PPP(l) 

K=I ABS (M-M) 

K= (N-M  )/K 

I-H 

M I N  =  M I 

NIN*NI 

D=S (N)-S(M) 

IF  (D*(SI (NI  )-SI  (MI  )  )  )  10,20.20 

MIN*NI 

NIN=MI 

KI=  IABS(NlN-MIIv) 

IF  (KI)  <tO»^0,30 

KI«  (MN-MIN)  /K  I 

II=M1N-KI 

C  =  0. 

IF  (MODE)  60*60,50 

C  =  l. 
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60 


7  0 


6  0 


IOC 


1 1  =  1 1  +  K I 

S  S  =  S  I  (  1  I  ) 

I=I  +  K 

IF  (1-K)  tfO*93*8G 

IF  (0*(i>(  I)-SS  )  )  70*70*  90 

J  =  I 

I-I-K 

SS=  SS-S  ( I  ) 

FPPPP=C*(FPPP(J)-FPPP(I))/(5(J)-S(I)) 

F  F  =  F  P  P  P  (  l)  +  .25*S5*Ff'PFP 

FF=FPP (I  )+  iS*FF/2. 

FF=FP(A)+.i*SS*FF 

FI(  ID  =F(  I)  +SS*FF 

IF  (II-NIN)  60*100*60 

CONTINUE 

RETURN 

END 


1C 


subrgu 

lkJO*  XS 

GENERA 

COMMON 

1*A0(12 

2  *  C 1  ( 1 7 

32(17)* 

<»HA»  C  A* 

DIMENS 

N-l 

IF  (XS 
DZ=Z(K 
IF  (PA 
IF  (PS 
SCALP= 
IF  (PS 
TX=3.0 
SX=-SC 
IF  (  IP 
CALL  P 
CONTIN 
IPLDT= 
CALL  F 
CALL  P 
ENCODE 
CALL  S 
ENCODE 


TINE  THREED  ( I  PLOT*  S  V*  SM*  CP*  X  »  Y*  T  IT  LE*  OC  »  A  L*  ZONE*F  M2,Q  EV  *  C  HO 

CAL* PSCAL ) 

TEj     THPtE     DIMENSIONAL     CALCOMP     PLOT;>     ON     CnC     6t00 

G{129*2  6*17)*SC(129»17)*E0(17)*IV(129,17)*ITE1(17).ITE2(17) 
*)*Ai(129)*A2(129)*A3ll29)*BC(26)#31(26),62(26)*B3(26)*Z(17) 
)*C2(17)*C3(17)*CM17),C5(17)»XC(17)*XZ(17)»XZZ(17)»YC(17)*Y 
rZZ(17)»KSYM,NX*r>if,NZ*KTfcl,KTE2,l5YM,SCAL,SCALZ*XK,GivEGA,ALP 
SA*F1ACH*T»PA*ABLA0E 
ION     X(l)»     Y(l),     S\M1)*     SM(D*     CF(1),     TITLE  (10),     R(20) 

CAL.Nt .0.)     SCALX».5*A8S(XSCAL)/CH0RD0*ZONE 

TE2)-Z(KTE1 ) 

.tQ.l. )     DZ=1./C5(KTE2)-1./C5(KTE1) 

CAL.GE.U.)     SCAL>=5./DZ 

-1.00 

CAL.ME.O.)     SCALP=-.5/ABS(PSCAL) 

AL>*XC  UTE1  ) 

LOT. ME. 1)     GO    TO    10 

LOTSdL     (  1000*25hANT0INE     BOURGEADF    +  337wv»H) 

UE 

C 

R  AME 

LOT     (1.25*1. »-3) 

(65, 190, R)     FMACH*F;i2*DEV*AL 
YMBOL     (0.0*0.75*  .  1  <• »  R*0.*65) 

(bO*200*R) 
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CAlL    SYrttlQL     (  .  50,  1 .25,  .  14, R,  0.  ,60) 

CONT INUE 
20  CONTINUE 

K  =  l 
30  CONTINUE 

K  =  K  +  1 

IF  (  K  .  G  T  .  •<  F  E  2  )  GO  TO  70 

IF  (K.LT.KTE1)  GO  TO  3  3 

U-ITEl(K) 

I2-ITE2U) 

CALL  V  E  L G  (K,K,SV»$IVI>CP,X,Y) 

SY=5.*(Z(K)-Z(Krfcl))/(Z(<TE2)-?(KTtl))+2.45 

SCP  =  5.*(Z(K)-Z(KTFl))/<Z(«rt2)-Z(KT::l))+2.75 

DO  4C  1=  11 » 12 

X ( I  ) =iC ALX*X  ( I  )  +  SX 

Y( I  ) =  SC  ALX*Y ( I  )  +SY 

CP(  I)=SCALP*CP  (  I  )  +SCP 
4  0  CON! I  Nut 

IF  (M.EJ.2)  GO  Tu  50 

N=I2-I 1+1 

CALL  LINt  (<(ll)»CP(Il).N,i,C,2,0.,l.,0.»l.) 

GO  TO  30 
50  CONTINUE 

N= I  2-1 1  +  1 

DC  60  1  =  11 f  12 

X ( 1  >  =  X  (  I  )  +  TX 
60  CONTINUE 

CALL  LINf 

GO  TU  2C 
70  CONTINUE 

M  ■  M  + 1 

IF  (M.GT 

GO  TO  20 
80  CONTINUE 

CALL  FRAME 

CALL  VELC 

Il=ITEi(KTEl  ) 

12= I TE2 ( KTE 1  ) 

SCALX=1.5*XK/3 

DO  90  I'll. 12 

X(  1  )  =  2.*X(  I) +SCALX4SX  +  1. 

Yd  )=2.*Y(  I  )*SC  ALX+  .5 

CP (1  )  =  Y(  I  )+6. 
90  CONTINUE 

N=I2-I1+1 

CALL  LINE  (  X(  I  I  )»  Y(  I  I  )  ,  N  ,  1  ,  U  »  2  ,  0  . ,  1  .i 

CALL  LINE  (X{I1)iCP{I1)»Nj1iOi2»0.iL 

ENCODE  (60. 210, R) 

CALL  SYMBOL  (1.0, 3. 5. .14, P.O. ,60) 

ENCODE  (60.220, »)  DC 

CALL  SYMBOL  ( 1 . 0, 7 . 75 , . 1 4, R , 0 . , 60 ) 

CALL  PLOT  (-1 . 25,-1. »-3) 


(X(ll),f(Il),N,l,0,2,0.,l.,C.,l.) 


2)  GO  TO  80 


(KTE1,KTE1,SV»SM,CP»X,Y) 


14159263 


1  .  ) 
.1.  ) 
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1CC 


LX=NX/2+l 

K  Y  =  N  Y  +  1 

N  A  =  P  A 

KD=  («Tt2-ME  1+1  )/2*NA  +  (  l-NA)*(KTfc2-<TEl+2  ) 

KD=l"iAX0(  K.D  i2) 

K  =  K  Ttl 

C0N1 I N  U  c 

CALL    FRAKE 

CALL     PLOT     (  3.  ,  <• .  J>,-3) 

DO    110     J  = 2  * K  Y 

0=30(J)+SQ(LX»k) 

A   _    v   /-     /    ^     *     •      *    I      Tr     /     i     [I    C 


G=30(J)+bQ(LX»K) 

A=XC(K)  +  AL3G( ABS ( 1 .-. t*SCAL*U*D) )  /XK 

IF     (0*LJ.LT.T)     GO     TO     12C 
110    CONTINUE 
120    CONTINUE 

0=B0(KY)+SD(2»K) 

B  =  l.  +  ( Au(2)*A0  12  )-0*0  )/T 

j=SCal*Ao(2)*j 

B=XC(K)+ALJG(3*B+0*0)/XK/2. 

A=-2./ A 

B  =  <t.  /B 

IF      (A.Lt.C.)     A  =  6 

SSX=1.2/Z0NE*XK 

SSX=AMIN1(SSX»A»6) 

DO  130  I=2>NX 

A=AO(  I  ) 

LPL0T=3 

DO  13  0  J=?»KY 

0=BO(J)  +  SO(1M) 

X1*1.+.5*SCAL* (A*A-0*Q) 

Y1=SCAL*A*G 

X2=XC(*<)+AL0G(Xl*xl  +  ri«ri)/xK/2. 

Y2=YC (K )+ATAN2 (Y1>X1 ) /XK 

Y3=bIN( Y* )/C5 (K  ) 

X2=SiX*X^ 

Y3=SSX*Y3 

CALL  PLOT  (X2,Y3>LPL3T) 

LPLCT  =  2 
130    CONTINUE 

LPLQT=3 

DO     140    J=2»KY 

0=80( J)+SO(LX, K ) 

XI  -  1 .-.5*SCAL*D*0 

X2*XC(K)+ALQG(ABS(X1)  )  /XK 

Y2  =  YC ( K )-ATAN2 (0. »  XI) /XK 

Y3»SIN( Y2)/C5(K  ) 

X2=SSX*X2 

Y3=SS**Y3 


Y3=SSX*Y3 

CALL    PLOT     (X2#Y3*LPL0T) 

LPL0T=2 


l^C    CONTINUE 

DO     160    J=2,KY 
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LPLCT=i 

CO    loG    1=2 >N< 

A=AC( I ) 

D=dO(J)+SO(I»<) 

X1  =  1.  +  .5*SCAL*  (  A* A-J*0  ) 

Y1=SCAL*A*H 

X2  =  X C  (  K  )  +A LUG  (  X 1  *  V  1  +  f  1  *Y  i  )  / XK  /  2  . 

Y2=YC(KJ+ATAN2(Y1»X1)/XK 

Y3=SIN(Y2)/C5(K) 

X2=S$X*X2 

Y3=SSX*Y3 

IF     (  i . N  E . L  X )     3 C     10    15  0 

Y2=2.*YC(K)-Y2 

Y^=SSX*SIM( Y2 ) /C5 (K ) 

CALL    PLOT     (X2>Y4»LPL3T) 

LPLUT=i 
150  CONTINUE 

CALL  PLOT  (  X2»  Y3»  LPLuT ) 

LPL0T=2 
16C  CONTINUE 

Zl*PA*ZJNt*EX?(Z(<))+(i.-DA)*Z(K) 

ENCODE  (bO,230,R)  Zl 

CALL  SfMbLL  (  -  1.  0.  4.  5  ,  .  14,  »,  C  .  ,  60  ) 

CALL  PluT  (-3.  ,-<t  .  5,-3  ) 

K=K+kO 

IF  (K.LE.KTE2)  GO  TO  100 

K  =  K-1 

IF  (K.FC.KTE2)  GO  TO  100 

DOl./XK 

IF  (PA.NE.1.0)  GO  TO  160 

N  =  2 

CO  170  1=1, N 

S  =  FLOAT(  1) /FLJAT(N) 

X0=0.5/Z0NE*S 

CALL  CUT  (XO,OC,SCALPj5V,SM,CP,X,Y,Z,YC»ITE1,1TE2,KTE1,<TE2,KSYMJ 
170  CONTINUE 
160  CONTINUE 

RETURN 


19  0  FOR 
16HA 
2  00  FOR 
210  FDR 
220  FOP 
230  FOR 
END 


CM     {  5HM1     =     t  F4  .2.  1H»  2*  >  5HM2    =     ,  F  5  .  2  t  1  H,  2  x  ,  oH  DE  V     =     ,F5.1,H,2X> 

LP    =     ,F*.l  ) 

MAT  (21HPR6SSURE  DISTR  I  BUT  ION  *  5  X  ,  l^HAL  A3E  prfjFlLt  ) 

MAT  (<:3H  CASCADE  REPRESENTATION) 

MAT  (9H     G/C  =,F5.2) 

MAT  (25H  GRID  ON  THE  iUPFACt  I     =»F5.2) 


10 


jG 


*0 


5C 


GO 


s 

IK 
1 
C 

l1* 

D 

i  Y 

c 

N 
[ 
C 
I 
I 
t: 
D 
1 
C 
1 
E 
C 
F 
S 
C 

D 
M 
I 
C 
B 
C 
F 

D 
R 
P 

C 

c 
c 
c 

C 

D 
x 
D 
x 
Y 
C 
C 
D 
A 
X 


UBkQuTINE  CLT  (XC,DCj;>CALP,SV>:>''i>CP>X,Y»Z,YC,ITFl,ITE2,l<Ttl,'<TE2j 

SYM 

HIS  SlrJKOJTlNE  PLOTS  5ECTIl>n3  OF  The  CD'IPi*  t  S  SO" 

CMMJN  /S1''/  jl(129»2fc)»G2(l29»26)*SX(129)»$Z(129),':>XX(129)>iX2U2 

>»SZZ(129),SJ(129),*1(129)>C(129)>D(1?9)  ,  1 1  ,  i  P  ,  L  X  .  'i  x  ,  \  y  ,  v.  Y  ,  1  L  ,  A  A  C 

01,  J2,f^ 

I M  E  N  S  1 0  M  S  V  (  1  )  »     S  M  (  1  )  »  C  P  ( 1 )  >     X(l),  Yd),  HI),     I  T  t  i  (  i  )  ,     I  T  t  2  (  1  )  , 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Department,  nor   any 
person  acting  on  behalf  of  the  Department: 

A.  Makes  any  warranty  or  representation, 
expressed  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information   contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or   process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus 
method,  or  process  disclosed  in  this 
report. 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Department"  includes  any  employee  or 
contractor  of  the  Department,  or  employee  of 
such  contractor,  to  the  extent  that  such 
employee  or  contractor  of  the  Department,  or 
employee  of  such  contractor  prepares,  dissemi- 
nates, or  provides  access  to,  any  information 
pursuant  to  his  employment  or  contract  with  the 
Department,  or  his  employment  with  such 
contractor. 
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